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1 .  INTRODUCTION 

This  report  summarizes  a  three  year  effort  under  AFOSR  contract 
to  elucidate  plausible  velocity  coupling  mechanisms  in  solid 
propellant  rocket  instability. 

"Velocity  coupling"  is  a  rather  general  term,  which  would  tend  to 
encompass  all  phenomena  not  included  in  pressure-coupled 
instability?  in  practice,  the  majority  of  instability  phenomena 
in  solid  propellant  rocket  motors  seem  to  resist  explanation  by 
pure  pressure-coupled  instability.  In  most  general  terms, 
pressure  is  a  state  variable,  and  the  response  of  the  propellant 
combustion  to  pressure  perturbations  could  be  explained  without 
detailed  account  of  the  associated  fluid  dynamic  processes  within 
the  motor  —  provided  the  gaseous  combustion  region  were 
sufficiently  thin  to  be  conceived  as  quasi-steady  (in  the  sense 
of  responding  "instantaneously"  to  any  pressure  perturbation). 
This  leaves  the  thermal  wave  relaxation  within  the  condensed 
phase  as  the  only  possible  mechanism  for  dynamic  coupling  with  an 
external  acxoustic  field,  considering  a  Rayleigh  instability 
mechanism.  This  is  in  general  insuf fucient.  Fluid  dynamic 
phenomena  must  be  incorporated  to  account  for  other  instability 
mechamisms.  This  means,  in  particular  that  the  instability 
characteristics  of  a  given  propellant/motor  configuration  can  not 
be  conceived  globally,  based  on  the  propellant  properties  alone: 
one  *ust  take  into  account  the  particular  internal  flow  field,  or 
fluid-dynamic  variables.  This  complicates  the  task  of 
instability  calculation  considerably,  for  obvious  reasons. 

After  the  effect  of  pressure  perturbation  has  been  separated, 
what  seems  to  remain  is  the  (external,  tangential)  velocity 
effect;  hence  the  name  for  the  instability  coupling  mechanism. 

In  practice,  this  may  involve  turbulence-combustion  interaction, 
residual  exothermicity  (slow  kinetics)  coupling,  visco-acoustic 
coupling,  vortex-boundary  interaction,  or  other  diverse 
phenomena,  all  having  to  do  with  instability  mechanisms  coupled 
with  the  fluid  dynamics  of  the  internal  flow. 

The  phenomena  under  consideration  are  also  characteristically 
nonlinear,  in  the  sense  that  the  response  depends  upon  the 
perturbation  amplitude,  and  efficient  coupling  (and  energy 
exchange)  exists  between  different  frequency  components  (in 
particular,  high-frequencies  influencing  low-frequency  behavior). 
This  accounts  for  the  considerable  complication  in  the  analysis, 
and  renders  limited  utility  to  quasi-linearized  calculations. 

Against  this  general  background,  the  present  investigation  has 
engaged  in  the  following  three  major  tasks,  undertaken  during  a 
period  of  three  years:  (1)  a  critical  literature  evaluation 

regarding  mechanisms  of  velocity-coupled  instability,  (2)  a 
detailed  order  of  magnitude  evaluation  of  the  various  mechanisms, 
and  (3)  a  detailed  numerical  analysis  of  the  most  prominent 
mechanism. 
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The  major  accomplishments  of  the  present  study  are  in  the  three 
aforementioned  areas:  The  critical  literature  review  included 
vortex/boundary  interaction,  an  area  studied  previously  by 
Flandro,  Combustion/Acoustic  interaction  in  the  coreflow,  due  to 
residual  exothermicity ,  studied  previously  by  Ben-Reuven  and 
Caveny  in  conjuction  with  nitramine  rocket  propellants,  and 
visco-acoustic  coupling.  All  of  these  areas  have  been  covered 
extensively  in  the  two  foregoing  annual  interim  scientific 
reports  under  this  contract. 

The  subject  of  turbulence-combustion  interaction  has  also 
received  due  coverage  during  the  present  study,  although  outside 
of  the  areas  covered  in  the  two  previous  interim  reports. 

An  additional  objective  of  the  study  was  to  establish  and 
maintain  contact  with  the  experimental  study,  carried  out  in 
parallel  at  UTC/CSD,  under  Dr.  Robert  Brown.  This  has  been 
likewise  carried  out,  through  numerous  discussions  and  exchange 
of  information,  which  definitely  shaped  both  the  focus  and  the 
outcome  of  the  present  velocity  coupling  analysis. 

Aside  from  the  literature  surveys  published  earlier,  cf  AFOSR  TR- 
82-1017,  this  analysis  has  two  major  accomplishments:  (1) 
Identification  of  the  visco/acoustic  coupling  as  a  potentially 
powerful  velocity-coupled  instability  mechanism,  through  both 
detailed  order  of  magnitude  study,  as  well  as  small-perturbation 
(asymptotic)  analysis,  and  (2)  development  of  an  explicit 
numerical  finite-difference  algorithm,  for  solution  of  the 
nonsteady,  compressible,  axisymmetric  flowfield  within  the  rocket 
motor,  geared  toward  the  elucidation  of  the  visco-acoustic  effect 
predominant  near  the  injected  surface. 

An  interesting  outcome  of  the  present  study  is  in  the  area  of 
steady  injected  internal  flow.  The  singular  small-perturbation 
analysis  has  enabled  analytical  expressions  for  the  dimensionless 
axial  pressure  drop,  and  the  surface  friction  coefficient,  in 
terms  of  the  injection  Reynolds  number  and  Mach  number.  These 
afford  exceptionally  good  agreement  with  the  experimental  data 
obtained  by  Dr.  Brown  at  CSD/UTC  recentlu,  as  well  as  earlier 
measurements  by  Olson  and  Eckert.  These  tend  to  be  collaborated 
by  the  numerical  results  generated  recently  by  this  author.  The 
evidence  points  to  the  plausibility  of  the  visco/acoustic  driving 
mechanism  (as  a  necessary  condition);  in  the  meantime,  one  need 
not  invoke  turbulence  or  other  phenomena  to  explain  the  axial 
pressure  drop  deviation  from  the  pure  inviscid  (rotational) 
theory  of  Culick  Taylor,  Berman,  and  others.  The  reader  is 
referred  to  Appendices  A  and  B  for  detailed  discussions. 

The  small  perturbation  boundary  layer  analysis  of  the  viscous 
sublayer  at  steady  state  has  been  the  subject  of  two 
publications,  presented  at  the  recent  (1983)  20th  JANNAF 
Combustion  Meeting,  and  at  the  Aerospace  Sciences  Meeting  of  the 
AIAA  (Jan.  1984);  the  latter  has  been  submitted  to  the  AIAA 
Journal  for  publication.  Both  papers  are  enclosed  in  Appendices 
A  and  B  mentioned  above. 
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At  the  closing  of  the  present  research  period,  the  computational 
finite-difference  algorithm  developed  has  just  begun  to  yield 
useful  results.  Numerical  stability  was  demonstrated  in  marching 
toward  steady  state  with  fixed  injection  rate  and  wall 
temperature  imposed;  the  initial  data  utilized  Culick's  inviscid 
solutions.  Some  of  the  results  following  1200  timewise 
integration  steps  are  shown  in  Figures  2.1  through  2.4  herein. 

In  particular,  the  computed  axial  pressure  drop  shown  in  Figure 
2.1,  is  similar  to  that  measured  by  Brown  et  al  (1983),  and 
correlated  by  the  perturbation  analysis  of  Ben-Reuven,  all  of 
which  depart  from  the  axial  pressure  drop  predicted  by  the 
inviscid  theory.  Integration  can  be  readily  carried  out  with 
imposed  time-dependent  boundary  data,  to  simulate  the  perturbed 
exit  nozzle  behavior  in  the  CSD  experiment. 

The  explicit  finite  difference  algorithm  developed  herein  follows 
a  modular  design,  and  allows  for  arbitrary  radial  mesh 
specification  (within  the  limits  of  physical  resolution  and 
numerical  stability).  The  modular  structure  can  readily 
accommodate  additions  [more  comprehensive  boundary  data 
treatment;  accounting  for  combustion;  including  turbulence 
transport] ,  in  the  form  of  subprograms.  A  joint  publication  with 
Prof.  Vichnevetsky  is  in  preparation  on  the  numerical  stability 
features  when  source  terms  are  incorporated. 

The  remainder  of  the  present  report  comprises  a  User's  Manual  for 
the  "MOSCO"  finite  difference  program.  Three  chapters  provide 
background  information  regarding  the  formulation  of  the  partial 
differential  system  and  the  finite  difference  algorithm  used  for 
solution,  a  discussion  of  the  program  features,  and  a  discussion 
of  the  results,  respectively. 
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2.  FUNCTIONAL  DESCRIPTION  OF  THE  FINITE  DIFFERENCE  CODE. 

2.  1  INTRODUCTION. 

MOSCO  is  a  finite  difference  FORTRAN  code,  to  solve  the  viscous, 
time-dependent,  compressible  Navier-Stokes  equations.  The  four 
equations  of  motion,  for  an  axi -symmetr ic  geometry  are: 
continuity,  axial  momentum,  radial  momentum  and  energy.  The 
finite  difference  algorithm  used  is  a  modified  unsplit  MacCormack 
scheme*  The  general  mode  of  integration  is  explicit  predictor- 
corrector. 

For  a  complete  description  of  the  original  scheme  and  a 
discussion  of  the  limitations  thereof,  the  reader  is  referred  to 
the  two  original  MacCormack  papers  listed  in  References  1  and  2. 

A  detailed  discussion  of  the  partial -differential  system  is  given 
in  Ref.  3  herein.  It  should  be  pointed  out  that  for  a  nonsteady, 
compressible  configuration,  tha  flowfiald  ia  non-aolanoidal , 
i.e.,  divergence  of  the  velocity  vector  is  nonzero.  This  gives 
rise  to  a  number  of  terms  emanating  from  the  symmetrical  part  of 
the  stress  tensor,  which  can  not  be  ignored  in  general. 

The  two  important  modifications  made  to  the  original  scheme  are 
as  follows:  (u  only  advactive  terms  are  written  in  the 

conservation  form  (general  flux  terms),  while  dissipative  terms 
of  all  kinds  appear  as  "source"  terms  ( non-conser vat l ve )  in  the 
finite-  differenced  equations  of  motion.  Also,  <ii>  arbitrarily- 
varying  radial  mesh  size  is  employed,  along  with  a  uniform  axial 
mesh;  along  with  the  variation  imposed  in  (1),  this  neccess itates 
proper  representation  for  the  second-order  and  radially-space- 
centered  first-order  derivatives,  to  reflect  the  aforementioned 
nonuniformity  and  maintain  second-order  consistency  of  the  scheme. 

Actual  implementation  of  the  MOSCO  program  involves  very  small 
radial  mesh-size  near  the  r=l  boundary,  representing  the  porous, 
injected  wall.  The  mesh  is  then  gradually  increased  toward  r=U, 
at  the  cylindrical  centerline.  A  ratio  of  1/100  between  the 
respective  mesh  sizes  has  been  tried.  A  similar  ratio  would  also 
hold  normally  between  the  smallest  radial  increment  and  the 
(uniform)  axial  mesh  size. 

Preliminary  small  perturbation  analysis  has  pointed  out  the 
necessity  of  proper  accounting  for  the  near-wall  processes,  where 
viscous  dissipation  is  dominant.  The  importance  of  these  viscous 
processes  to  proper  understanding  of  the  steady  state  axial 
pressure  distribution  has  been  discussed  elsewhere.  The 
interaction  between  acoustic  vibrations  and  the  dissipative 
processes  in  this  nearfield  drives  the  visco-acoustic  coupling 
mechanism,  an  important  component  of  nonlinear  velocity  coupled 
instability.  For  this  purpose,  the  fine  radial  stepsize  near  the 
wall  is  implemented. 


The  unsplit  MacCormack  scheme. 


under  these  conditions. 


has  one 
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major  drawback :  to  maintain  numerical  stability,  the  Courant- 
Fr 1 edr 1 chs- Lewy  condition  is  implemented  to  calculate  the  time 
step  sice  according  to  the  smallest  spatial  increment. 

Cosequent 1 y ,  the  overall  (uniform)  tima  *tap  is  very  small, 
considering  the  value  implied  by  the  radial  mesh  near  the 
centerline,  or  the  regular  axial  mesh.  Over  a  period  of  several 
thousands  of  timesteps,  this  might  lead  to  an  appreciable  error 
regarding  the  timewise  profiles.  In  this  respect,  the  split 

MacCormack  method  might  be  superior  for  obtaining  long-term 
transient  or  low  frequency  behavior. 

The  present  algorithm  is  set  up  to  integrate  the  fully  dynamic 
system.  Its  limitations  regarding  long-time  and  low  frequency 
behavior  are  recognized.  The  major  advantage  of  the  system 
relative  to  implicit  integration  algorithms  is  quite  obvious,  in 
obtaining  typically  low-cost  and  low  CPU  performance.  Other 
advantages,  due  to  the  specific  constr uct ion ,  are  the  modular 
structure,  in  which  changes  in  the  formulation  are  easy  to 
implement,  (details  specific  to  the  partial  differential  system 
being  solved  appear  only  in  three  subprograms)  additions  can  be 
readily  made  (through  auxilliary  subprograms)  and  programming 
errors  can  be  detected  easily. 

The  general  structure  of  0SC03,  the  cold-flow  simulation  FORTRAN 
code  is  depicted  in  Fig.  1.  This  is  a  development  based  on  the 
earlier  ROSCO - ser les ,  reported  previously  by  this  author. 

The  logic  of  the  timewise  integration  cycle  is  shown  in  Fig.  2. 
The  particular  configuration  corresponds  to  the  marching  toward 
steady  state  mode,  pertaining  to  the  configuration  of  the  program 
listing  in  Appendix  A  herein. 

A  concise  description  of  the  partial  differential  system,  the 
boundary  data,  and  the  finite  difference  method  is  given  in  the 
following  Section.  This  is  followed  by  a  detailed  glossary  of 
input  parameters,  and  the  output  maps.  Options  regarding  the  use 
of  output  data  files  were  not  included;  other  options  which 
require  extensive  setup  (and  depend  crucially  on  printer 
configuration),  such  as  printer-plots,  were  omitted.  The  program 
listing,  and  a  sample  output  are  given  in  Appendices  A  and  B, 
r  espect i ve 1 y . 
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The  question  of  relevance  of  cold-flow  simulation  to  internal 
rocket  flows,  and  the  use  of  porous  materials  to  simulate 
propellants  has  been  raised  many  times.  Based  on  the  present 
experience,  it  seems  that  both  are  highly  relevant,  although 
obviously  defficient  in  several  ways.  For  instance,  normal 
propellants  respond  to  a  positive  pressure  perturbation  by 
raising  the  burn  rate,  while  a  porous  plug  would  yield  a  lower 
mass  flow  rate  if  the  pressure  differential  across  it  decreased. 
The  inability  to  simulate  the  high  thermal-gradient  region  near 
the  surface  of  a  propellant  with  the  existing  injected  flow 
techniques  is  another  drawback.  On  the  other  hand,  probing  within 
the  interior  of  a  burning  propel] ant  grain  is  still  quite 
impossible.  The  conclusion  drawn  is  that  cold  flow  injection  and 
porous  materials  have  many  limitations,  but  could  still  be 
successfully  utilized  for  study  of  certain  limited  physical 
interactions,  certainly  not  the  entire  instability  phenomena.  An 
excellent  example  of  such  limited  application  is  the  visco¬ 
acoustic  mechanism,  which  comprises  a  clear  fluid  dynamic 
phenomenon . 

The  importance  of  interior  fluid  dynamic  considerations  to  the 
proper  understanding  of  solid  propellant  motor  behavior  can  not 
be  overemphasized.  Whether  it  is  the  highly-ordered,  frequency 
dependent  visco-acoustic  phenomenon,  acoustic-combustion 
interaction,  or  similar  mechanisms  of  interaction,  or 
turbulence/combustion  effects--  all  require  details  of  the  local 
fluid  dynamic  processes.  It  is  therefore  strongly  recommended 
that  both  theoretical  analysis  and  idealized  experimental 
simulation  of  the  internal  flow  field  of  solid  propellant  rockets 
be  continued. 


2.6  CONCLUSIONS  AND  RECOMMENDATIONS 


A  considerable  amount  of  experimental  data  has  been  generated  by 
the  cold-flow  simulation  effort  at  UTC/CSD,  by  Dr.  Brown. 
Nonsteady  wall  heat  transfer  coefficient  measurements  as  well  as 
hot  wire  anemometry  were  utilized.  Since  the  wall  heat  transfer 
is  very  important  to  the  understanding  of  combustion  instability, 
the  folowing  comments  are  made.  It  appears  that  the  coherently 
averaged  spectra  of  the  wall  heat  transfer,  which  filter  out  the 
RMS  or  DC-effects  due  to  turbulence  and  viscous  second-order 
interactions,  demonstrate  the  validity  of  the  visco-acoustic 
coupling  mechanism,  suggested  by  this  author  in  our  1982/S3 
Interim  Report.  The  original  figures  from  Brown's  Annual  Report 
have  been  merely  re-grouped  here,  in  Figs.  2.8  through  2.12. 

These  figures  show  that  the  effect  of  increasing  peturbation 
amplitude,  (from  A' /A  of  0.6  to  5.9*)  and  increasing  perturbation 
frequency  (from  76  to  170  Hz)  both  cause  enhanced  wall  heat 
transfer,  coherently  with  the  acoustic  frequency,  at  axial 
stations  considerably  upstream  of  the  velocity  antinode  of  the 
first  or  second  axial  acoustic  modes.  All  of  these  effects  are 
explained  by  the  visco-acoustic  coupling  mechanism. 


When  the  same  data  were  re-cast  with  the  DC  effects  included,  it 
is  evident  that  the  low  amplitude  phenomena  are  swamped  by  near¬ 
wall  turbulence  effects  and  other  high  frequency  noise,  peaking 
at  the  surface  (as  corresponding  hot-wire  traverses  indicate). 
This  occurs  at  roughly  10  diameters  downstream  of  the  head  end. 
The  higher  perturbation  aaplitude  and  frequency  coherence, 
nevertheless,  persist  way  downstream  of  this  point. 

Other  hot  wire  anemometry  data,  obtained  by  Dr.  Brown  during 
1983,  seemed  to  indicate  relatively  high  turbulence  intensities 
near  the  centerline  (close  to  the  head  end)  while  these 
intensities  increased  toward  the  wall  at  downstream  stations.  It 
should  be  emphasized  that  these  data  were  mostly  obtained  at  the 
low  perturbation  amplitude  of  A' /A  =  0.6*  ;  close  to  the 

centerline,  the  axial  velocities  are  high,  but  du/dr  is 
vanishing;  in  the  meantime,  both  the  radial  velocity  component  as 
well  as  dv/dx  are  vanishingly  small.  Near  the  wall,  however,  the 
radial  velocity  is  nonzero,  and  both  du/dy  and  dv/dr  are 
appreciable,  while  the  axial  velocity  is  small.  The  foregoing 
physical  picture  holds  for  most  of  the  flowfield.  One  therefore 
expects  the  vorticity  effects,  responsible  for  turbulence 
generation,  to  be  strongest  near  the  walls,  even  though  the  walls 
are  injected. 

It  is  therefore  suggested  that  the  results  pertaining  to  this 
low-amplitude  of  perturbation  be  thoroughly  checked  prior  to  any 
far  reaching  conclusions  regarding  the  origin  or  evolution  of 
turbulence  in  internal  injected  flows.  The  experimental  evidence 
of  most  significance  (and  relevance  to  velocity  coupling 
instability)  seems  to  be  the  steady  state  data,  and  the  high 
perturbation  amplitude  and  frequency  data. 
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On  the  other  hand,  the  corresponding  (not  normalized)  radial 
velocity  distribution  at  the  same  station,  shown  in  Fig.  2.7, 
departs  appreciably  from  its  inviscid  counterpart.  This  is  not 
surprising,  as  it  infers  that  du/dx  (axial  acceleration)  is 
somewhat  less  than  the  inviscid  prediction,  since 

~  +■  v/r ") 

according  to  the  continuity  equation  with  negligible 
compressibility;  so  that  whenever  dv/dr  obtains  a  positive 
increment,  as  occurs  near  the  surface  in  Fig.  2.7,  du/dx 
decreases  accordingly.  This  is  consistent  with  the  information 
obtained  in  the  axial  pressure  drop.  Fig.  2.4,  which  shows  that 
the  computed  p-drop  is  somewhat  smaller  than  the  theoretical 
inviscid  prediction,  since  the  axial  acceleration  balances  the 
axial  pressure  gradient  in  the  coreflow,  cf  Appendix  A. 


In  conclusion,  the  steady  state  results  herein  are  both  self- 
consistent  as  well  as  in  agreement  with  the  predictions  of  the 
viscous  ( perturbational )  wall  layer  analysis,  and  the  recent 
measurements  of  Brown.  The  algorithm  has  been  demonstrated 
numerically  stable  in  this  mode,  and  can  be  utilized  for 
nonsteady  flowfield  simulation. 
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2.5  DISCUSSION  OF  RESULTS 

Some  of  the  results  obtained  from  the  numerical  solutions  will  be 
discussed  herein,  pertaining  to  marching  toward  steady  state.  The 
reader  is  referred  to  Figs.  2.4  -  2.7  in  which  converged  data 
following  1201  timewise  integration  steps  is  summarized. 
Considerably  more  data  is  available  at  both  the  last  step  and 
intermediate  time-levels  as  well;  a  partial  output  listing  of  the 
same  datum  case  discussed  herein  is  appended  to  the  User's  Guide, 
Chapter  3  . 

It  should  be  pointed  out  that  there  is  nothing  particular  about 
the  Datum  Case  configuration  used.  The  specific  radial  mesh 
divisions  are  in  no  way  "the  best"  or  optimal.  The  axial  mesh 
size  is  quite  large  (1.0  in  the  dimensionless  system,  i.e.,  equal 
to  one  radial  length);  smaller  mesh  size  and  considerably  more 
points  in  both  radial  and  axial  directions  would  definitely  yield 
better  quality  of  data.  The  computer  program  can  accommodate  both 
with  only  minor  (and  obvious)  modification  to  Common  statements 
and  the  input  data. 

The  physical  input  data  used  was  to  simulate  cold  nitrogen 
injection  under  experimental  conditions  similar  to  those  used  by 
Dr.  Brown  at  CSD,  withou  attempting  to  simulate  any  particular 
set  of  conditions  exactly.  Again,  conditions  like  the  injection 
velocity,  temperature,  chamber  reference  pressure,  etc  can  be 
readily  varied,  over  a  sufficiently  wide  parametric  range. 

In  the  meantime,  it  is  asserted  that  the  finite  difference  scheme 
is  numerically  stable,  within  the  limits  of  time-resolution 
discussed  in  Section  2.3,  and  that  the  formulation  as  well  as  its 
implementation  correspond  to  the  compressible  Navier-Stokes 
equations,  in  axisymmetric  form.  In  numerous  tests  with  coarse 
radial  mesh,  (the  viscous  scales  being  much  smaller),  numerical 
instability  consistently  developed;  the  aforementioned  stability 
therefore  obviously  depends  on  the  degree  of  spatial  resolution 
near  the  surface. 

In  Fig.  2.4  he  axial  pressure  drop  is  plotted  against  axial 
distance,  x.  the  trend  in  departure  from  the  inviscid  solution  of 
Culick  is  the  same  as  that  measured  by  Brown  at  CSD.  In  Fig.  2.5, 
the  friction  coefficient  and  Stanton  number  (heat  transfer 
coefficien)  are  plotted  against  the  injection  ratio.  The  nearly- 
linear  trend  is  similar  to  that  measured  (for  Cf)  by  Olson  and 
Eckert.  Departure  from  the  inviacid  theoretical  normalized  axial 
velocity  profile  is  quite  small  at  10  radial  distances  from  the 
head  end,  as  shown  in  Fig.  2.6. 
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2.4  THE  BOUNDARY  CONDITIONS 


As  differs  from  the  physical  plane,  all  of  the  dependent 
variables  in  the  numerical  simulation  must  be  specified  along 
each  of  the  four  boundaries.  The  available  physical  boundary  data 
has  been  discussed  already  in  section  1.2  herein. 


Along  the  centerline,  (x,r=0>,  all  of  the  boundary  data  is  physically 
available.  The  radial  gradients  of  the  density,  axial  mass  flux, 
and  enthalpy,  are  all  zero;  the  radial  velocity  is  zero: 


U(l,  1,  K)=U(1,£,  K) 
U<3,  1 ,  K)  =U  < 3,  £,  K) 
U  <4,  1,K)=U<4,£,K) 


>  d(23i,K)  =  0. 


Along  the  porous  sidewall,  (x,r=l)  the  axial  velocity  is  zero, 
the  injection  velocity  and  the  wall  temperature  are  specified, 
and  the  pressure  extrapolated,  using  a  three-point  algorithm: 

C - POROUS  SURFACE:  PRESSURE  EXTRAPOLATED.  WALL  ENTHALP=HWP. 

DO  3  K=£, KXX 

U<4,  JRR,  K)  =  (RRS*U<4, JRM1,K)-U(4, JRMl-i, K> ) /RSM1 
3  U(l,  JRR,K)=U<4,  JRR,K)/HWP 

Es*»1^ae.s-1  . 

This  allows  calculation  of  the  density  at  the  wall: 


P  ~  f°/hu)  (t^oop  #3  A WVZ)  . 


At  the  head  end,  <x=0,r>,  both  velocity  components  are  zero, 
while  the  wall  temperature  is  specified.  The  pressure  is 
extrapolated  axially  to  the  wall,  using  a  two-point,  second  order 
accurate  extrapolation: 

C - HEAD  END : X=Q.  NOTE:U£=0  AND  U3=0,  ALWAYS.  WALL  TEMP. =HW=F < R > . 

DO  £  J*1,JRM1 

U<4,  J,  1>=U<4,  J,  £>*£.  — U  (4,  J,  3) 

U(l,  J,  1>=U(4,  J,  1)/HW0 
DO  £  M*l,  MXX 

2  U(M,  J,  KXX)  =U<M,  J, KXM1 )*£.  -U (M,  J,  KXM1 - 1 > 

From  which  the  density  at  the  wall  is  calculated: 

k*=<OAi«  xrC^Tji) 

At  the  exit  plane,  all  of  the  variables  are  extrapolated  axially, 
using  a  two-point  formula: 

00,1000=  2-UttKXX-O-  0 (T, KXX-2-) - 

which  simulates  a  continued  duct,  not  the  entrance  to  a  choked 
nozzle;  to  simulate  the  latter,  characteristics  segments  have  to 
be  utilized  locally,  in  a  finite  difference  form. 
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cur  rent  1  y  used,  and  the  (regular)  axial  mesh  size  are  such,  that 
the  minimal  timestep  is  always  obtained  by  use  of  delta-R, 
namely : 


At*  =  aAv/Vwc (c -tu)  < Ate. 


(20) 
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'djQ  _i  JL_  (  _  Qr^-i  C 

2rr±  AG  ■(&}-,  L  A ^  Afj-i  )k. 


(  18) 


Obviously,  the  differences  between  adiacent.  radial  increments 
should  be  small  where  the  third  order  radial  derivative  is 
expected  to  be  large,  to  insure  small  truncation  error  for  the 
second  derivative.  This  principle  is  indeed  followed  in  the 
neighborhood  of  the  porous  wall;  near  the  centerline,  however, 
the  third-order  radial  derivatives  are  expected  to  be  quite 
negligible,  so  that  a  larger  radial  meshsize  variation  can  be 
implemented . 


The  reason  for  the  peculiar  weighting  factors  in  the  radial  flux 
terms  within  the  predictor  and  corrector  steps  is  now  clear,  to 
maintain  second-order  accuracy  in  the  overall  (combined)  timestep 
integration;  otherwise,  (if  the  same  finite-differencing  is  used 
for  both  axial  and  radial  derivatives),  the  algorithm  is 
consistent  only  to  first  order,  with  a  radial  diffusion  term 
which  is  proportional  to  the  difference  between  two  adjacent 
radial  increments. 


The  timestep  size  is  determined  by  the  smallest  spatial 
increment,  using  the  Courant-Friedrichs-Lewy  stability  condition , 

At5AV  =  ,  c=Vmc  ,19, 

JjIC 

where  the  coefficient  B  <  1  (strictly),  is  typically  taken 
between  O.b  and  0.6. 

Previous  error  analysis  of  a  modified  Lax-Wendroff  scheme  (one 
dimensional,  nonste  dy),  which  involved  source  terms  has 
indicated  that  using  a  CFL  number  B=0.b  has  clear  advantages 
toward  reduction  of  the  induced  numerical  diffusivity.  Due  to  the 
inherent  similarities  between  both  mathematical  systems 
(predominantly  hyperbolic  in  axial  direction),  and  numerical 
method  used  (overall  central -di 1 ferences  in  space,  predictor- 
corrector,  explicit),  the  same  CFL  number  is  currently 
implemented.  Performance  in  marching  toward  steady  state  is 
clearly  superior  to  identical  cases  using  datum  case  input  with 
B-0.8,  in  terms  of  improved  stability.  In  this  mode  of 
operation,  the  cost  is  an  obviously  longer  overall  itegration 
time. 

The  relative  magnitudes  of  the  smallest  racial  increment 


2.3  THE  FINITE  DIFFERENCE  ALGORITHM:  2-STEP/UNSPLIT  MacCORMACK 


Following  the  original  scheme  by  MacCormack ,< 1971 ,  1975)  the 
following  explicit  two-step  procedure  is  proposed,  in  a 
predictor-corrector  manner: 


A£i _ 

^0  4AVJ-l 


ta  !• 


At  S, 


vt 


(13) 


xr^ 


J.  at. 

2  v  a«~- 


( 14) 


where  overbar  denotes  "predicted”  properties.  The  mam 
differences  from  the  original  MacCormack  scheme  are  in  the  source 
terms  herein  (the  original  was  written  for  conservation-form 
differential  system),  and  the  special  treatment  required  by  the 
uneven  mesh.  The  "source"  terms  herein  contain  the  dissipative 
effects  in  the  system. 


Fig.  3  depicts  the  spatial  mesh  in  the  axi -symmetric  flowfield. 
The  configuration  employs  a  nonuniform  radial  mesh  distribution, 
which  is  very  fine  near  the  porus  wall  (r=l>  and  much  coarsec 
near  the  centerline  (r=0>.  The  axial  mesh  is  uniform. 


The  overall  configuration  of  the  two  steps,  (first  backward,  then 
forward)  when  combined,  is  space -centered .  To  remain  at  the 
second-order  acuracy  level  in  the  radial  direction,  for  the 
source  or  S-terms  which  involve  first  and  second  derivatives,  the 
following  finite  differencing  algorithm  is  utilized: 

^ Q  /Sen- G^-iy1  £Wi 
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transonie  region  is  outside  of  the  domain  of  interest  for  the 
present  system. 

The  initial  data  is  specified  as  the  following  arbitrary  (in 
principle)  distributions, 

V(t =0,  r,x)  =V°CrJ  aO  U2) 

The  foregoing  formulation  is  the  subject  for  the  numerical 
OSCG/COLD-FLOW  algorithm.  Obviously,  mathematical  closure  is 
obtained  when  additional  boundary  data  are  invoked:  this  is 
attained  by  using  the  appropriate  characteristics  relations 
near  the  boundaries.  The  auxiliary  boundary  data  are  discussed  in 
the  following  section. 


The  radial  and  axial  flux  terms  are,  respectively, 


ft  =  5>v2?  f u-Vj  "Yf hv) 

0T  *  (?“,  tu2+wc-  ,  rfhu  3 

The  source  terms  are  defined  as  follows.  Si  =  0, 

O  -  J/'gV  J_^F  . 

°z  r  V-sr  p  /  2nr  “r 

f  •* +!P]-  ho 


+  cg^f} 


(8) 


The  following  physical  boundary  data  are  available  for  the  cold- 
flow  simulation  problem;  on  the  centerline,  (t,  r=0,  x) 


\!^0 


9 


ur 


=  o 

3  dr  s> 


(9) 


at  the  injected  porous  surface,  (t,  r=l,  x) 


v/--VoCxJt)^  vu-0  J  k  =  hu/CxJt^)  (10) 

at  the  nonpermeable ,  solid  head-end,  <t,  r,  x=0> 

V-G  7  \k~0  ?  k  =  hH(rjt}  (id 

The  exit  plane,  defined  by  (t,  r,  x=L),  forms  an  entrance  into  a 
short  convergent  nozzle  section,  assumed  choked  at  all  times: 
detailed  description  of  the  nozzle  entrance  and  its  inherent 
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2.2  FORMULATION  OF  THE  DIFFERENTIAL  SYSTEM 


The  four  equations  of  motion  pertain  to  compressible,  viscous, 
nonsteady  motions  in  an  axisymmetric  flowfield.  The  system  is 
written  in  differential  form,  employing  source  terms  to  represent 
the  various  dissipative  effects. 

The  following  dimensionless  independent  variables  are  introduced, 

r=r*/eS,  x  =  xVe,*,  t  - tVt;  a> 

where  to  =  Rq/vq. 

The  corresponding  dependent  variables  are: 

V=V'/vo'e  ,  u-  =  U*/VoX 

f  =<f  /?„*  ,  k  =  KV hi  ,  =fVP.‘r  <2> 

The  properties  used  for  non-dimensional ization  are  the  injection 
velocity,  the  reference  gas  density,  and  the  reference  chamber 
pressure.  The  corresponding  thermal  enthalpy,  hG  ,  is  calculated 
from  the  caloric  equation  of  state, 

where  gamma  =  Cp/Cv  is  considered  a  constant.  The  reference 
adiabatic  speed  of  sound  is, 

<  =  (7  fz=  \Jc^ovT 

The  corresponding  injection  Mach  number,  Reynolds  and  Prandtl 

numbers  are,  respectively, 

M.  e 

^  ,  Vr^ljfcf/X 

Thus,  the  dimensionless  differential  system  can  be  written  in 
general,  for  0  <  x<L,  0  <  r  <  1,  and  t  >  0  : 


and  the  dependent  variable  vector  is: 

C  T ,  fv  ,  p=fK) 


ChxIO^  SURFACE  STANTON  NUMBER 


USER ' S  MANUAL  for  "MOSCO" 
COLD  F~LOW  SIMULATION  PROGRAM 


'  v  tmm*  «4*V*-  ,A*C*A  _ 


3.1  INPUT  DESCRIPTION 

3.1.1  BLOCK  DATA  Statement. 


BLOCK  DATA 

COWON/AREAO/VZERO,  PSTAR,  TSTAft,  VISC,  COND,  UBAR,  RUJ,  HSTAR,  CP, 

>CFl,  RSTAR,  XSTAR,  XO,  IPRNT1,  IPRDT  (2),  IC0RR(2) ,  INC  (3),  ITHAX,  VO,  P0.M 
COtNON/AREAl  /GAMA,  6AH1, GAH2,  GAH3,  GAH3T,  C81,  CS2,  CG3,  C64,  C65, 

♦REO,  EHO,  EH02,  PI,  PIT,  EKH,  RH06,  SSND,  PRN,  DX,  TOOX,  DX2,  RR<26) ,  DRRI25) , 
*R2 (26) , R3 (26) ,  XX (26) ,  DT,  EPS 

COWOM/AREA2/HXX,  JRR,  KXX,  JRM.KXM,  KXPR(6) ,  KXPtfX,  JRPR  (6), 
♦JRPMAX,  DGX (4, 26, 26) ,  DFR(4, 26, 26) , VIC<4, 26, 26) , HUO, HUP, RRS, RSM1 
DATA  COW),  PSTAR,  RSTAR,  XSTAR,  GABA,  CFL/.  8166, 2.  NES, .  85, .  55, 1. 4, .  5/ 
DATA  VZERC, VISC, TSTAR, UBAR, RUJ/1.000,  l.E-5,  278.,  .828,  8.3l4/ 
DATA  JRPBAX,  JRPR,  JRR/6,  1,3,5,7,9,11,  12/ 

DATA  KXPMAX,  KXPR,  KXX/6,  1,3,5,7,9,11,  12/ 

DATA  VO, PO.RO/1.0000, 1.0000,1./  , ITHAX, IPRNT1/1201, 200/ 

DATA  IPRDT,  1C0RR,  1MC/-1,0,  0,1,  -1,0,1/,  HXX/4/ 

DATA  DRR/. 05,.  15,. 2,  .15, .15,  .l,.l,  .06,. 025,. 01,. 005,  14«0./ 

END 


All  of  the  input  enters  through  the  BLOCK  DATA  statement,  and 
transmitted  to  the  remainder  of  the  program  through  the  labelled 
Common  Blocks. 

Throughout  the  program,  all  dimensional  parameters  are  in  SI 
units.  The  integration  itself  is  carried  out  with  dimensionless 
variables;  suitable  transformation  of  the  differential  system  has 
been  made,  prior  to  utilization  of  the  finite  difference 
algorithm  herein. 
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The  seven 
glossary . 

Line  1 

COND 

PSTAR 

RSTAR 

XSTAR 

GAMA 

CFL 

Line  2 

VZERO 

vise 

TSTAR 

WBAR 

RUJ 

Line  3 

JRPMAX 

JRPR 

JRP 

Line  4 

KXPMAX 

KXPR 

KXX 


input  data  lines  are  discussed  in  the  following 


thermal  conductivity,  J/m-s-R  of  the  gas  (air) 
reference  pressure,  N/m^ 
inner  radius,  m 
axial  length,  m 

ratio  of  specific  heats,  CA/CV 
Courant-Fr iedr icks-Lewy  number 


injection  velocity,  m/s ec 
viscosity  coefficient,  Kg/m-s 
ref.  temperature  of  gas,  K 
mean  molecular  weight  of  gas,  kg/mol 
universal  gas  constant,  J/mol-R 


=  total  number  of  points  to  print  out  in  radial 
mapping  (see  output) 

=  radial  position  vector;  numbers  correspond  to 
distinct  radial  stations,  and  their  total  number 
should  agree  with  JRPMAX. 

=  maximum  number  of  meshpoints  in  radial  direction 
(including  ends) 


(same  as  JRPMAX)  total  number  of  axial  points  in 
printout  of  axial  maps 

axial  position  vector;  axial  mesh  indices  of 
points  to  be  printed  out;  total  number  should  be 
equal  to  KXPMAX 

maximum  number  of  axial  mesh  points,  including 
boundaries 
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Line  5 


VO 

= 

dimensionless  (absolute)  injection  velocity 

PO 

= 

dimensionless  reference  pressure 

RO 

= 

dimensionless  inner  channel  radius 

ITMAX 

= 

total  number  of  timewise  integration  steps 
allowed 

IPRNT1 

= 

number  of  timewise  integration  steps  between 
output 

Line 

_6 

IPRDT 

= 

predictor  step,  index  displacement  vector 
(backward  differences),  used  in  DGDX  and  DFDR 

ICORR 

= 

same  as  above,  for  corrector  step  (forward 
differences) 

INC 

= 

index  shift  vector  for  the  central-difference 
algorithms  of  SORCE  subprogram 

Line 

_7 

DRR 

— 

radial  increment,  starting  with  increment  near 
centerline,  and  ending  with  that  near  the  wall 
(last) 

NOTE:  DRR  is  dimensioned  to  26  and  empty  places 

above  12  must  be  padded  with  zeros. 

CAUTION:  Performance  of  the  program  is  quite 
sensitive  to  (small,  required)  step  size  near 
wall . 

This  concludes  the  entire  input  data  set. 


n  n  n 
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3.2  MAIN 


- MBR  VERSION  OF  NACCOfWCK  INTEGRATION  SCHEME.  12/9/83 

- FIRSr  TRIAL:UNSPLIT,  SINGLE  PREDICTQR/CQRRECTDR  CYCLE  EACH  BT. 

- VARIABLE  RABIAl  STEPSJZE.  ARBITRARY  DR(J)=INPUT  VECTOR. 

COWWAREAO/ VZERO,  PSTAR,  TSTAR,  VISC,  COND,  K8AR,  RUJ,  HSTAR,  CP, 

+CFL,  RSTAR,  XSTAR,  XO,  iPRNTl, IPRDTI2), 1C0RRI2), INC (3), ITNAX,  VL,  PG,  Ru 
COWION/AREA1/GANA,  SAMI,  GAME,  GAH3, 6AM3T,  CS1,  C62,  C63,  CS4,  CSS, 

♦REQ, EMO,  EHG2,  PI , PIT, EKM,  RH06, SSND,  PRN, CX,TODX,  DX2,  RR126),  DRR 125) , 
+R2 (26) , R3I26) ,  XX (26) ,  DT,  EPS 

COWWN/AREA2/MXX,  JRR,  KXX,  JRN1 ,  KXMl, KXPR (6) ,  KXPflAX,  JRPR  (6) , 
♦JRPflAX,  DGX  (A,  26, 26),  DFR  (4, 26, 26),  VIC  (4, 26, 26),  HUO,  HUP',  SRS,  RSH1 
C0MM0N/AREA3/LU4, 12, 12),UE(4, 12, 12),S(4, 12, 12),SB(4, 12, 12) 

TIM£=0. 

CSL.  SCAT  A 
K* 

DO  1  L=1,ITMAX 
IF(H.LT.  IPRNTDGQ  TO  2 
CALL  PRINT2(TI*E,L> 

K=« 

2  COfflNUE 

CAu  TIMNT1  (MXX, JRR, KXX, Tittt, DT) 

H=*+l 

1  COV1NUE 
S'0P2222 
END 


The  MAIN  program  just  calls  for  calculation  of  the  fixed 
parameters  and  the  initial  data  profiles,  all  in  subroutine 
SDATA.  Then,  in  Do-loop  #1  timewise  integration  is  called  out 
each  timestep,  by  calling  subroutine  TIMINT.  Printing  of  the 
variable  tables  of  interest  Is  made  by  calls  to  sobroutine  PRINT, 
each  preselected  number  of  timesteps:  the  integer  IPRINT1 
controls  the  number  of  timesteps  between  such  output  dumps. 
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3^3 _ TIME-INTEGRATION  PROCEDURE 

3.3.1  Subroutine  TIMINT 


SUBROUTINE  TIMNT1 IMF, JF,KF, TINE, DT> 

C - M6R  VERSION  Of  UNSPLIT  NACCORMAIK  TINE  ARCHING.  12/9/63. 

COHMCN/AREAQ/VZERO, PSTAR,  TSTAR,  VISC,  COND,  WBflR,  RUJ,  HSTAR,  CP, 

+CFL,  RSTAR,  XSTAR,  XO, 1PRNT 1, IPRDT (2) , ICORR (2) , INC (3) , ITMfiX, VO, PO,  RO 
COWON/AREA2/MXX,  JRR,  KXX,  JRMl ,  KXK: ,  KXPR  (6) ,  KXPMflX,  JRPR  (6 ) , 
♦JRPMAX,  D6X  (4, 26, 26) ,  DFR  (4, 2£,  26) ,  VIC  (4, 26, 26) ,  HUO,  HwP,  RRS,  RSM1 
COWI3N/AREA3/UU,  12, 12),UB<4, 12, 12),  S  A, 12, 12>,SB<4, 12, 12) 

C - PREDICTOR  STEP 

CALL  DFDR  (MF,  JF,  KF,  U,  2,  IPRDT) 

CAlL  DGDX (MF, JF, KF, U, 2, IPRDT) 

CALL  SORCEKMF,  JF,KF,U,S) 

DO  6  M=1,MF 
DC  6  J=£, JRM1 
DO  6  K=2,KX*tl 

6  L‘R ( •?,  J,  K)  =Li («,  J, K/-DT* (I>-'R («,  J, K)  +D6X  («,  Jf KJ -S <K,  J, X) ) 

C - CORRECTOR  STEP 

Cft-L  EAT)RY<Mf,Jf,KF,UB) 

CA. L  IFDR (!*!-,  JF, KF, UB, 2, ICORR) 

BL.  DGDX (MF, J^,Kr,UB,2, ICORR/ 

CAw_  SORCE 1  INF,  JF,  KF,  UB,  SB ) 

DO  7  *=1,MF 
DO  7  J=2,  JRK1 
DO  7  K=2, KXM1 

U(K,  J,K)  =  (U(M,  J,K)+UB(M,  J,K))/2.-.5*D7*(DFR(K,  J,K)+DGX(N,  J,K) 
-«S(H,  J,K)+SB(M,  J,K)  )/2. ) 

7  CONTINUE 

Cft..  BNDRY  (KF,  JF,  KF,  U) 

TIME=TIME>DT 

RETURN 

END 


Subroutine  TIMINT  carries  out  the  timewise  integration  process 
for  the  unsplit  MacCormack  scheme.  The  Predictor  step  involves 
calculation  of  the  radial  and  axial  flux  terms,  based  on  backward- 
differences:  this  is  carried  out  for  the  entire  interior  of  the 
physical  field  by  calls  to  subprograms  DFDR  and  DGDX, 
respectively.  The  integer  vector  IPRDT=(-1,0>  controls  the 
backward  differencing.  A  call  to  subroutine  SDRCE  establishes  the 
"source"  terms,  involving  dissipative  effects.  The  second  order 
differential  terms  are  discretized  as  central  differences  in 
subroutine  SORCE.  Subsequent  to  these  calls,  the  entire 
dependent-var iable  (predicted)  vector,  UB,  is  calculated  in  Do- 
1 oop  #6 . 


Continued 
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RHOG=PSTftR*WBftR/ (RUJ*TSTAR) 

SSN3=  (GAWkPSTAR/RHOG)  ##.  5 

REO=RfOG*VZERD*RSTAR/VISC 

EPS=l./SQ«T(REO) 

E*3=VZER0/SSKD 

E*l32=EM0w3 

XC=XSTflR/RSTAR 

GA*!]=6AtfH. 

6AX2=6AK1/GANA 

GA.y2=6AMfi*EBG2 

6Pi'3T=6M3/2. 

CS:=l.  /SfW3 
262=4. /3.  /RED 
263=1. /RED 
CG5=SAfll*6AH3/RE0 
**=:./(GAM3«R£0) 
HSTflR=PSTAR/<GAK2#RHOG) 
CP=HSTAR/TSTAR 
PRN=VISC*CP/COND 
C64=G(W/RE0/PRN 
JRfl=JRR-l 
KJ«i=Kn-i 
CM  ;=+5XX-l 
DX=XG/KXW1 
T0DX=2.*DX 
DX£-DX*DX 


Continued 


The  first 
RHOG 
SSND 
REO 

EPS 

EMO 

XO 

EKM 

HSTAR 

CP 

PRN 

GAM*  » 


group  of  parameters  calculated  are  as  follows, 
ref.  gas  density,  kg/m3 

ref.  adiabatic  speed  of  sound  in  gas,  m/s ec 

injection  Reynolds  number,  based  on  the  injection 
velocity  and  the  chamber  radius. 

the  small  parameter  pertaining  to  viscous  effects 
injection  Mach  number 
dimensionless  axial  length 

second  dimensionless  parameter  denoting  the  ratio  of 
inertial  to  viscous  effects:  Km  in  the  analysis  of 
Appendix  A. 

ref.  specific  thermal  enthalpy,  J/kg 
ref.  specific  heat,  isochoric,  J/kg-K 
Prandtl  number 

constant  parameters  associated  with  gamma,  the  ratio  of 
specific  heats. 


CG  *  • 


constant  parameters  associated  with  Reynolds  number 


-46- 


3.5  SUBROUTINE  SDATA  -  CALCULATION  OF  CONSTANT  PARAMETERS  AND 
THE  INITIAL  DATA 


SUBROUTINE  SDATA 
DIMENSION  VR(26),VX<26) 

LOGICAL  DIFFPR 

COMMON/ AREAO/VZERO,  PSTAR, TSTAR, VI SC, COND, NBAS, RUJ, HSTAR,  CP, 

+Ca,  RSTAR,  XSTAR,  XO,  IPRNT1,  IPRDT  (2) ,  ICORR  (2),  INC<3>,  UMAX,  VO,  PO, » 

C0MM0N/AREA1/GAMA,  SAMI ,  GAME,  GAM3,  GAM3T, CS1, C62,  C63,  CS4,  CB5, 

+RE0,  EMO,  EM02,  PI, PIT, EKM, RHOG, SSND, PRN, OX,  TODX, DX2, Rfl (26) ,  DRR (251 , 

♦R2(26),R3(26),XX(26),DT,EPS 

C0MM0N/AREA2/MXX, JRR, KXX,  JRM1, KXM1, KXPR (6) , KXFW)*, JRPR (61 , 

+JRPI1AX,  D6X  (4, 2b,  26) ,  DFR  (4, 2b,  26) ,  VIC  (4, 26, 26),HWO,  HUP,  RRS,  RSM1 
COMMON/ AREA3/U  (4, 12, 12),UB(4, 12, 12),S(4, 12, 12),SB(4, 12, 12) 

C0MN0N/AREA4/D1FPR 

Continued 


All  of  the  information  calculated  herein  is  passed  to  the 
remainder  of  the  program  through  the  labelled  common  blocks. 


Continued 


3.4.2 _ Subroutines  XPRNT1  and  XPRNT2 


SUBROUTINE  XPRNT 1  HE, JF,  KF,  UP,  ZX,  ZY) 
DIMENSION  UP<NF,JF,26),ZX<JF),ZY<26> 

1H  F0RMAT(/,2X,'RH0  =  DIM.LESS  DENSITY') 

WRITE (6, 111) 

WRITE(6,102XZX<N),N=1,JF) 

DO  1  K=1,KF 

1  WRITE (6, 103)ZY(K),  (UP(1,  J,K),J=1,JF) 

lW  FORMAT (/,2X,* VR  *  DIN. LESS  RADIAL  VELOCITY') 
WRITE (6, 184) 

WRITEI6,102)(ZX(N),N=1,JF) 

DO  2  K=1,KF 

2  WRITE (6, 1*3)  ZY (K) ,  (UP (2,  J,  K) ,  J=l,  JF) 

1*5  FORMAT  </,2X,'VX  =  DIN. LESS  AXIAL  VELOCITY') 
WRITER  105) 

WRITE(6,102)  (ZX(N),N=1,JF) 

DO  3  K=1,KF 

3  WRITE (6, 103)ZY(K), (UP(3,  J,K), J=l,  JF) 

106  FORMAT (/,2X,’P  =  DIN.LESS  STATIC  PRESSURE’) 
WRITE  (6, 106) 

WRITE  (6, 102)  (ZX(N),N=1,  JF) 

DO  4  K=1,KF 

4  WRITE (6, 103)ZY(K),  (UP(4,  J,K),  J=l,  JF) 

102  FORMAT (6X, ' YZ' , 6(2X, ' XZ=' , F7. 4) ) 

103  F0RMAT(2X,F6.3,6(1PE12.4)) 

RETURN 

END 


SUBROUTINE  XPRNT2(ST,CF,FLX,XX,KXX) 

DIMENSION  ST (KXX) , CF (KXX), FLX (KXX), XX (KXX! 

WRITE (6, 201 ) 

201  FORMAT <//,2X,’ SKIN  FRICTION (CF)  t  STANTON  NO. (ST)’,/) 
WRITE<6,202) 

202  FORMAT <6X,  ’  X’ ,  8X,  ’  CF’ ,  10X,  ’  ST' , 2X, ’ U1 <W) /U2 (CL) ’,/ > 

DO  204  1=2, KXX 

WRITE(6,203)XX(I),CF(I),ST(I),FLX(1) 

204  CONTINUE 

203  FORMAT (2X,  F5. 2, 5(  1PE12. 4) ) 

RETURN 

END 
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C - NOTE:  M  VR  AND  DVR  ARRAYS  ARE  TRANSPOSED  (ROTATED  ABOUT  THE 

C - MAJOR  AXIS)  RELATIVE  TO  THE  PARENT  U-ARRAY.  THIS  IS  ONLY  FOR 

C - PRINTOUT  REASONS 

WRITE (6, 98) TIME, L 

C - X-MAK 

WRITE  (6, 100) 

CALL  XPRNT 1 (MXX, JRPHAX, KXX, VX, RP, XX) 

IF(DIFFPR)GO  TO  5000 
WRITE (6, 201 )L 
WRITE (6, 100) 

CALL  XPRNT  1 (MXX, JRPMAX, KXX, DVX, RP, XX) 

CALL  XPRNT2 (ST, CF, FLXI, XX, KXX) 

C - R-MAPS 

5000  WRITE  (6,200) 

CALL  XPRNT 1 (MXX,  KXPMAX,  JRR, VR,  XP,  RR) 

IF (DIFFPR)  GO  TO  5001 
WRITE(6,201)L 
WRITE (6, 200) 

CAL.  XPRNT!  (MXX, KXPMAX,  JRR, DVR,  XP,  RR) 

5001  CONTINUE 

100  FORMAT <  //,  2X,  ’ X-MAP,  AT  DISTINCT  RADIAl  POSITIONS’,/, IX, 75(’-'),/> 

200  FDRMAT(//,2X,,R-MAP,  AT  DISTINCT  AXIAL  POSITIONS’,/, IX, 75!'-’),/) 

201  FORMAT (//, 2X, ’ DIFFERENCES  BETWEEN  CURRENT  TIMESTEP  N0.',H,/,2X, 
*’f*D  THE  INITIAL  DATA  U(M, J,KHJINITIAL(M, J,K)  s'  ) 

90  FORMAT  (//,2X,’TIfE  ,  1PE12. A, 5X, ’ TIMESTEP  NO.  =’,I4,/) 

RETURN 

END 


For  ease  of  carrying  out  the  printing,  the  radial  maps  VR(J,K> 
and  DVR < J , K )  are  transposed  (rotated  about  the  major  axis);  this 
allows  use  of  exactly  the  same  printout  procedure  in  subroutine 
XPRNT 1 ,  for  both  axial  and  radial  variable  maps.  Following  the 
appropriate  Table  Header  printouts,  the  four  consecutive  calls  to 
subroutine  XPRNT1  affect  printout,  with  four  table  groups,  each 
containing  four  single  tables  (one  for  each  variable).  The  cal  to 
subroutine  XPRNT2  affects  printout  of  the  friction  coefficient 
and  Stanton  number  Table,  following  the  X-maps. 

Note  that  the  length  of  each  table  is  not  limited;  the  only  self- 
imposed  format  limitation  used  in  the  design  is  the  compatibility 
with  a  printer  page  width  of  80  columns.  Consequently,  there  is 
no  attempt  to  fit  printout  groups  into  any  particular  page 
length . 
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DO  2M0  J=l,  JRR 
DO  2000  K=1,KXPW)X 
KK=KXPR(K) 

(MM)(1,J,KK) 

VR(l,K,J)=RHO 
VR(2,K,  J)=U(2,  J,KK)/RHG 
VR  (3,  K,  J)  =U(3,  J,  KK)  /RHO 
VR(4,K,  J)=0(4,J,KK) 

XP(K)=XX(KK1 
DO  2100  H=1,WT 

2100  DVR  (N,  K,  J)  =VR  (N,  K,  J)  -VIC  (N,  J,  KK) 

IF  (KK.  61.1)80  TO  2200 
DVR(3, 1,J1^0. 

GO  TO  2080  Continued 

2200  DVR(3,K,  J)=DVR(3,K,  J)/(FI«XP(K)) 

2000  CONTINUE 


i 


A  similar  procedure  is  carried  out  within  Do-loop  #2000,  for  the 
radial  printout  map,  VR,  and  its  corresponding  difference-map, 
DVR.  Distinct  axial  positions  are  now  used,  as  specified  in  the 
integer  vector  KXPR,  again  with  an  imposed  maximum  of  6 
positions:  for  each  of  these,  the  radial  dependent  variable 
distributions  are  printed  out,  with  XP  denoting  the  corresponding 
normalized,  distinct  axial  stations.  Note  that  in  this  instance 
the  radial  pressure  difference  map  (not  normalized)  is  calculated 
instead  of  the  aforementioned  axial  pressure  drop. 


Continued 
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C - LOAD  l-m A  ARRAYS 

DO  1000  K=l, KXX 
DO  1000  J=l,  JRPMflX 
JJ=JRPR(J) 

RHCftJd,  JJ,K) 

VXU,  J,K)=RH0 
VX  (2,  J,  K)  =U  (2,  JJ,  K)  /RHO 
VX  (3,  J,  K)=U<3,  J  J,  K)  /RHO 
VXU,  J,K)=U(4,  JJ,K) 

RP1J)=RR(JJ) 

«XfN4XX-l 
DO  1100  M=1,MXH 

1100  DVX («, J,K)=VX (H, J,K)-VIC(H, JJ, K) 

DVX  (4,  J,  K) = ( VX  (4, 1, 1 ) -VX  (4,  J,  K) )  /  ( VX  <4,  J,  K)  »EPS) 

IF (K.  BT. 1)60  TO  1200 

DVX(3,J,1)=0. 

GO  TO  1000 

1200  DVX (3,  J,  K) =DVX (3,  J,  K)  /  (PI«XX  (K) ) 

1000  CONTINUE 


Continued 


In  Do-loop  #1000.  the  primitive  dependent  variable  array  VX  is 
calculated  (still  dimensionless),  namely,  the  axial  and  radial 
fluxes  are  replaced  by  the  respective  velocities.  The  purpose  is 
to  create  an  X-map,  in  which  the  full  axial  distributions  of  each 
variable  are  printed,  at  distinct,  preselected  radial  positions; 
up  to  6  such  distinct  radial  positions  are  available  with  the 
current  setup,  to  facilitate  printing  with  an  80-column  nominal 
page  width.  The  integer  vector  JRPR  stores  the  distinct  radial 
positions  to  be  printed;  RP  denotes  the  radial  position  values, 
normalized . 

In  the  nested  Do-loop  #1100  the  difference  array  DVX  (departure 
from  initial  data)  is  calculated,  for  the  density  and  the  two 
velocities;  the  fourth  variable  calculated  therein  is  the 
dimensionless  axisl  pressure  drop,  normalized  by  the  small 
parameter  epsilon  (EPS).  The  axial  velocity  difference  is 
normalized  by  PI»X,  to  factor  out  the  linear  growth  (due  to 
accummulated  injection)  along  the  axis. 
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3_.4 _ OUTPUT  SECTION 

3.4.1 _ Subroutine  PRINT 


SUBROUTINE  PRINT2<TINE4L) 

DIMENSION  VK(4,6,26),M(4,6,2b),XP(6),VR(4,6,2b),DVR(4,6126) 

+  ,  RP  (6) ,  ST (26) ,  CF (26) , FLXI (261 

COMMON/ AREAO/VZERO,  PSTftft,  TSTftR,  VISC,  COND,  UBfift,  RUJ,  HSTAR,  CP, 

+CFL,  RSTAR,  XSTAR,  XO,  IPRNT1,  IPRDT  (2) ,  ICORR  (2) ,  INC  (3) ,  ITNAX,  VO,  PO,  RO 
COMMON/ftREfll /6ANA, 6AM1 , SANE, GAM3, 6AM3T, CS1 ,  C62,  CS3,  CS4,  CSS, 

+RE0,  EHO,  EN02,  PI, PIT, EXM,  RHOS,  SSND,  PRN,  DX,  TODX,  DX2,  RR  (26) ,  DRR  (25) , 

+R2(26),R3(26),XX(26)fDT4EfiS 

C0NM0N/AREA2/MXX, JRR,  KXX, JRM1, KXM1, KXPR<61 , KXPWTX,  JRPR (6) , 

+JRPMRX,  DGX  (4, 26, 26) ,  DFR<4, 26, 26) ,  VICI4, 26, 26) ,  HUO,  HUP,  RRS,  RSMl 
CQHN0N/AREA3/U (4, 12, 12) , UB(4, 12, 12) , S (4, 12, 12) ,  SB (4, 12, 12) 

C0MH0N/AREA4/DIFPW 
LOGICAL  D1FFPR 

C - CALCULATION  OF  SKIN  FRICTION  I  STWTON  NUMBER 

DO  5M  K=2,KXX 
RH01=U(1, JRM1,K) 

WALLH=U  ( 4,  JRR,  K )  /U  ( 1 ,  JRR,  K) 

CEN7H=U(4,!,K)/U(1,1,K) 

CF(K)=VISC#U(3,JRM1,K)/(RH01*DRR(JRM1))/(.5*U(3,1,K)«2/U(1,1,K)) 

ST  (K)=-COND/CP»  (U  (4,  JRM1,  K)  /RHOl-NALLH)  /DRR  <  JRN1) 

IF(CENTri.EQ.HALLH)H3  TO  5M 

ST (K) =ST (K) / (U (3, 1, K) * (CENTH-MALLH) ) 

FlXI(K)=U(2,JRR,K)/U(3,  1,K) 

see  CONTINUE  Continued 


Subroutine  PRINT  loads  the  output  maps,  calls  for  their  actual 
printing  (through  subroutines  XPRNT1  and  XPRNT2),  and  prints  the 
header  for  each  output  table. 

Skin  friction  coefficient  <CF>  and  heat  transfer  coefficient  (ST, 
Stanton  number)  at  the  porous  injected  sidewalls  are  calculated 
first,  using  2-point  differencing.  The  resulting  CF  and  ST 
vectors  will  be  printed  out  through  subroutine  XPRNT2,  vs  the 
axial  injection  (mass  flux)  ratio,  FLXI,  also  calculated  herein. 


Continued 
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3.3.5  Subroutine  BNDRY 


SUBROUTINE  BNDRY  Of,  JF,KF,U) 

DIMENSION  U(NF,JF,KF) 

C0PM0N/AREA2/HXX,  JRR,  KXX,  JRM1, KXM1, KXPRI6) ,  KXPMAX,  JRPR  (6) , 
♦JRPMAX,  D6X  (4, 26, 26) ,  DFR  (4, 26,  E6) ,  VIC(4, 2S,  26) ,  HHO,  HMP,  RRS,  RSN1 

c - MBR-VERSION  OF  BOUNDARY  DATA  TREATMENT. 12/13/63 

C - CENTERLINE,  R=8.  NOTE:  U2=«  ALWAYS. 

C - DU1 /DR=DU3/DR=DU4/DR=6  TO  MAINTAIN  AXIAL  SYMMETRY. 

DO  1  K=2*KB 
U(1,1,K)=U(1,2,K) 

U(3,1,K)=U(3,2IK) 

1  U(4,1,K)=U(4A« 

C - POROUS  SURFACE:  PRESSURE  EXTRAPOLATED.  NALL  ENTMLP=H*. 

DO  3  K=2,KXX 

U(4,  JRR,  K) = (RRS*U  (4,  JRM1 ,  K) -U  (4,  JRMH ,  K) ) /RSM1 
3  U(1,JRR,K)=U(4,JRR,K)/HHP 

C - ^  £ND;x=0.  N0TE:U2=8  AND  U3=8,  ALWAYS.  NALL  TEMP.  =rt*=F ( R) . 

DO  2  J-l, JRM1 

U(4,  J,  1  )=U(4,  J,  2)  *2.  -U  (4,  J,  3) 

U  ( 1,  J,  1)=U(4,  J,1)/HN0 
DO  2  M=1,MXX 

U(M,  J,  KXX  >=U  (M,  J,  KXM1  )*2. HI  (M,  J,  KXM1-1 ) 

2  CONTINUE 
RETURN 
END 


I* 


Subroutine  BNDRY  is  called  following  the  predictor  and  the 
corrector  steps,  from  subroutine  TIMINT.  The  mode  of  operation 
is  fully  explained  in  the  previous  chapter  of  this  manual. 


3.3.4 _ Special  Functions 


FUNCTION  DDZ(A,B,D) 

DDZ=(B-A)/D 

RETURN 

END 

FUNCTION  DDR2IA,B,C,D1*» 

C - SECOND  DERIVATIVE  FINITE  DIF.  -NONUNI FOM  *9* 

DDR2=2.  *  1DDZ  (B,  C,  D21-DDZ  <  A,  B,  D1  > )  /  ID1+D2) 

RETURN 

END 

FUNCTION  CDX2(A,B,C,D2) 

C - SECOND  DERIVATIVE  FINITE  DIF.  —UNIFORM  MESh,  CENTRAL  DIFF. 

CDX2=(ft-2.*B+C)/D2 

RETURN 

END 

FUNCTION  DBAR(A,B,C,D1ID2) 

DBAR= (D2*DDZ (A,  B, D1 ) +D1*DDZ (B, C, D2> ) / ID1+D2) 

RETURN 

END 

FUNCTION  DDRX1  (tF,Q,Dl,D2,DXT) 

DIMENSION  Q(NF,NF) 

C - MIXED-UP  SECOND  DERIV. ,  UNIFORM  X,  NON-UNIFORM  R. 

DU1=DDZ(Q(1,1),Q(1,3),DXT) 

DU2=DDZ(Q(2, 1),  0(2,3), DXT) 

DU3=DDZ  (0(3, 1),Q(3,3),DXT) 

DDRX1=DBAR(DU1,DU2,DU3,D1,D2) 

RETURN 

END 


C - LIBRARY  Of  FINITE  DIFFERENCE  FUNCTIONS: 

DVDX=DDZ ( VT  (2, 1 ) ,  VT  (2, 3) ,  TODX ) 

|  DUDX=DDZ (UT  (2, 1 ) ,  UT  (2, 3) ,  TODX) 

DPDX=DDZ (U tt,  J,  K-l ) ,  U<4,  J, K+l ) , TODX) 

DVDR=DBAR  ( VT  ( 1, 2) ,  VT  (2, 2) ,  VT  (3, 2) ,  DR1,  DR2) 
DUDR=DBAR(UT(1,2),UT(2,2),UT(3,2),DR1,DR2) 

DHDR=DBAR (HT ( 1 , 2) , HT <2, 2) ,  HT (3, 2) ,  DR1,  DR2) 

DPDR=DBAR  ( U  (4,  J-l,  K) ,  U  <4,  J,  K) ,  U  (4,  J+l,  K) ,  DR1,  DR2) 

|  D2VDX=CDX2 (VT (2, 1 ) , VT (2, 2) , VT<2, 3) , DX2) 

D2UDX=CDX2(UT(£,  1 ) ,  UT (2, 2) , UT (2, 3) , DX2) 

D2HDX=CDX2  (HT  (2, 1 ) ,  HT  (2, 2) ,  HT  (2, 3) ,  DX2) 

D2VDR=DDR2  ( VT  U ,  2) ,  VT  (2, 2) ,  VT  <3, 2) ,  DR1,  DR2) 

D2UDR=DDR2 (UT ( 1 , 2) , UT (2, 2) , UT (3, 2) , DR1 , DR2 ) 

D2HDR=DDR2(HT ( 1, 2) , HT (2, 2) , HT (3, 2) ,  DR1,  DR2) 

D2VDRX=DDRX1  (3,  VT,  DRl,  DR2,  TODX) 

D2UDRX=DDRX1  (3,UT,DR1,DR2,T00XI 
VOR=VT (£,  £)/«R(J) 

C - SOURCE  TERMS: 

S(l,J,K)=fl. 

S(2,J,K)=-CG1*DPDR  +  C62*(DVDR-V0R)/RR(J) 

»  +  +C63MD2VDX+D2UDRX/3. )+C62»D2VDR 

S (3,  J,  K)=C63* ( (DUDR+DVDX/3. >/RR(J) +D2UDR+D2VDRX/3. ) +C62f D2UDX 
S (4, J,  K) =SAM1* ( VT (2, 2) »DPDR+UT (2, 2) »DPDX) 

+  +C64* (D2HDR+D2HDX+DHDR/RR( J l ) 

♦  +CS5*  (2.  #  (DUDX»#2+DVDR«2*-V0R«2)  +  ( DUDR+DVDX )  **2 

-  (2.  /3. )  *  (DVDR+VOR+DUDX)  «2) 

|  |  •  i  CONTINUE 

RETURN 

END 


I 


Following  the  "LIBRARY  OF  FINITE  DIFFERENCE  FUNCTIONS",  first 
order  and  second  order  derivatives  are  calculated,  using  central 
differences  with  the  utility  functions  specially  constructed  for 
this  purpose.  The  variable  names  used  are  similar  to  the  actual 
differential  expressions.  This  form  allows  for  easy  detection  of 
errors . 

The  source  terms,  S(M,J,K>,  are  finally  calculated.  The 
parameters  involving  gamma,  injection  Mach  No,  and  Reynolds 
number,  are  all  pre-calculated  in  SDATA. 
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3.3.3  Subroutine  SORCE 


SUBROUTINE  SORCEl  <>F,JF,KF,U,9 
DIMENSION  U(MF,JF,KF),S(HF,JF,KF> 

DIMENSION  VT(313>lUT(3,3)lHTI3ia 

C - MBR  VERSION  OF  SOURCE  (DISSIPATIVE)  TERMS  IN  KCCORNACK  SDO£ 


C 

c 


COMMON/AREAO/VZERO,  PSTAR,  TSTAR,  VISC,  COND, H6AR, RUJ, HSTAR, CP, 

♦CFL, RSTAR,  XSTAR, XO, IPRNT1 , 1PRDTI2) , ICORR (2) , INC (3) , ITNAX,  VO,  PO,  RO 
COWON/AREA1  /GAMA,  SAM  1 ,  GAM2,  BAM3,  GAM3T,  C81 ,  CS2,  CS3,  CS4,  C65, 

+R£C’ DC'  ^  PI '  PIT>  SSNO,  DX,  TODX,  Die,  Rfi  (2b) ,  DfiS i 25), 

♦R2(26) ,  R3 (26) ,  XX (26) , DT, EJ5  ’ 

C0MH0N/AREA2/MXX, JRR, KXX, JRM1, KXM1, KXPR(6) ,  KXPMAX, JRPR (6), 
+JRPMAX,  D6X  (4, 26, 26) ,  DFR  (4, 26, 26) ,  VIC  (4, 26, 26) ,  HHO,  HUP,  RRS,  RSMl 
INC ( 1 ) =— 1,  INC(2)=«,  INC(3)=1  INDEX  SHIFT  OPERATI*. 

M=T)C  X-SHIFT  DUN. ,  N=T*  R-SHIFT  DOR. 

DO  I  K=2,KXM1 
DO  1  J=2,JRN1 
DO  3  M=l,3 
KK=K+INC(M) 


DO  3  N=l,3 
JJ=J+INC(N) 


RriO=U(I,  JJ,KK) 

VT  (N,M)=U(2,  JJ,KK)/RHO 
UT  (N,M)=U(3,  JJ,KK)/RH0 
3  HI  (N,M)=U(4,  JJ,KK)/RHO 
DRi=DRR(J-l ) 

DR2=DRR(J) 


Continued 


Subroutine  SORCE  carries  out  the  calculation  of  the  source  terms 
for  the  differential  system,  dissipative  terms  of  first  and 
second  order  are  calculated  by  use  of  central  differencing.  A 
dummy  index  shift  integer  vector  is  used,  INC= < - 1 , O , ♦ 1 )  ,  to 
facilitate  keeping  merely  an  array  of  3x3  of  each  dependent 
variable  (temporarily,  to  facilitate  the  local  calculations), 
these  are  the  local  arrays  RHO,  VT ,  UT  and  HT ,  with  (JJ.KK) 
serving  as  the  3x3  dummy  indices.  This  facilitates  great  savings 
in  core  and  storage  requirement. 


Continued 


SUBROUTINE  DFDR (MF, JF, KF, U, IF, JTR) 

C - HBR  VERSION  OF  RADIAL  FLUX  TERN  BENERATOR.  12/9/83 

DIMENSION  U(NF,JF,KF),JTR(IF),F(4,2) 

C0MM0N/AREA2/MXX,  JRR, KXX, JRM1.KXM1, KXPR<6>, KXPMAX, JRPR (6), 
+JRPMAX,  DGX  (4, 26, 26) ,  DFR  ( 4, 26,  £6) ,  VIC  (4, 26, 26) ,  HUO,  HUP,  RRS,  RSM1 
COMMON/ AREP1 /GAMA, 6AM1, GAM2, 6AM3, GAM3T, CGI, C62, C83, C84, CSS, 

♦RED, EHO,  EM02,  PI ,  PIT,  EKH, RH06, SSND,  PRN,  DX, TODX, DX2, RR (26) , DRR (25) , 
+R2(26),R3(26),XX(26),DT,EPS 
DO  1  K=2,  KXM1 
DO  1  J=2,  JRM1 
DRO=DRR(J+JTR(l) ) 

DR 1 =R3 ( J ) *DRO/DRR ( J- JTR  C  2 ) ) 

DO  2  Mil, IF 
JJ=J+JTR(N) 

F(1,N)=RR(J)*U(2,  JJ,K) 

V2=F(1,N)/U(1,  JJ,K) 

1F(N.GT.  l)GO  TO  6 

VM2=Vf 

GO  TO  7 

6  J1=JJ-JTR(U 

VM2=(RR(J1)*U(2,  J1,K)/U(1,J1,K)+RR(J)*U(2,  J,K)/U(1,  J,K))/2. 

7  F (2, N) =VH2*U (2, JJ,  K) 

F(3,N)=V2MJ(3,JJ,K) 

2  F (4, N) =V2«GANA*U(4, JJ, K) 

DO  3  H=1,MF 

3  DFR  (M,  J,  K)= (F  (M,  2)  -F  (M,  1 ) )  /DRl 
1  CONTINUE 

RETURN 

END 


Subroutine  DFDR  carries  out  the  radial  differencing  to  calculate 
the  radial  flux  array,  DFR.  Operation  is  very  similar  to  DGDX, 
with  JTR«<-1,0>  for  the  predictor  step,  and  JTR=<0,*1)  in  the 
corrector  step,  acting  as  the  appropriate  radial  index  shift.  The 
parameter  VM2  is  the  mean  radial  advective  velocity,  analogous  to 
the  UM2-term  discussed  earlier. 

The  differencing  is  not  symmetrical  overall,  as  the  parameter  DRl 
is  different  for  backward  and  for  forward  differences,  to  account 
for  the  variable  radial  mesh  size,  and  maintain  second-order 
accuracy  overall;  "overall"  herein  means  when  the  backward  and 
forward  terms  are  combined,  in  the  calculation  of  the  U-array 
following  the  predictor  step,  cf  Do-loop  #7  in  subroutine  TIMINT. 
When  a  uniform  radial  meshsize  is  imposed,  however,  the 
calculation  automatically  becomes  similar  to  the  overal 1 -central 
differencing  in  DGDX. 
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3.3.2  Subroutines  DGDX  and  DFDR 


SUBROUTINE  D6DX  Iff,  JF,  K F,  U,  IF,  KTR) 

C - MBR  VERSION  OF  AXIAL  FLUX  TERM  GENERATOR.  12/9/83. 

DIMENSION  U(HF,JF,KF),KTR(IF),G(4,2> 

COMMON/AREA1 /GAMA,  GAM1, GAM2, GAM3, GAM3T, 081,062,063,084,065, 

♦REO,  EHO,  EM02,  PI,  PIT,  EKM,  RH06,  SSND,  PRN,  DX,  TQDX,  DX2,  RR  (26) ,  DRR(25) , 
+R2 (26) , R3 (26) , XX (26) , DT, EPS 

COMM0N/AREA2/MXX,  JRR,  KXX,  JRM1.KXM1,  KXPR(6) ,  KXPMAX,  JRPR  (6) , 
+JRPMAX,  DGX  (4, 26, 26),  DFR  (4, 26, 26) ,  VIC(4, 26, 26) ,  HUG,  HUP,  RRS,  RSM1 
DO  1  J=2,JRM 
DO  1  K=2, KXM1 
DO  2  N=1,IF 
KK=K+KTR(N) 

U3=U(3,  J,  KK)/U(1,  J,KK) 

6(1,N)=U(3,  J,KK) 

G(2,N)=l)3*U(2,J,KK) 

IF  (N.  GT.  1)G0  TO  6 

UK3=U3 

GO  TO  7 

6  K1=KK-KTR(1) 

UM3=(U(3,J,Kl)/U(l,J,Kl)+U(3,J,K)/U(l,J,K))/2. 

7  G(3,N)=UH3*U(3,  J,  KK)  +CG1*U(4,  J,KK) 

2  G(4,N)=U3*SAMA*U(4,  J,KK) 

DO  3  H=l,ff 

3  D6X(K,J,K)=(6(M,2)-6(N,  I))/DX 
1  CONTINUE 

RETURN 

Ev'D 


Subroutine  DGDX  calculates  the  four  axial  flux  terms  in 
discretized  form.  KTR=<-1,0)  for  the  predictor  phase,  and 
KTR=(0,*1)  for  the  corrector  phase.  This  integer  vector  is  used 
as  an  index  transform,  to  calculate  the  dummy  axial  mesh  index 
KK.  This  way,  both  predictor  (backward  differencing)  and 
corrector  (forward  differencing)  can  be  facilitated  in  the  same 
subroutine . 

G(4,2)  is  the  temporary  advective  term  array,  reloaded  at  each 
new  meshpoint.  DGX  is  the  axial  flux  array  arising  from  this 
calculation . 

Note  that  the  use  of  the  Mean  axial  velocity  UM3  for  the  axial 
advective  term  G(3,n)  has  been  suggested  by  MacCormack ( 1974 ) ,  to 
alleviate  problems  inherent  in  the  non-conservative  nature  of  the 
RHO»U»U  terms. 


The  corrector  step  ia  similar  in  structure,  and  is  carried  out 
after  the  call  to  subroutine  BNDRY  establishes  the  variable 
boundary  data  corresponding  to  the  predicted  interior  UB- 
values.  Thus  calls  to  subroutines  DFDR,  DGDX  and  SORCE  establish 
the  corrected  flux  terms  and  source  terms  respectively.  The  value 
of  the  integer  vector  IC0RR=(0,1),  passed  in  the  argument  lists 
of  DFDR  and  DGDX,  insures  that  the  same  subroutines  will  carry 
out  forward  differencing  in  the  corrector  phase.  Do-loop  #7  is 
used  to  calculate  the  entire  inner  variable  array,  followed  by  a 
call  to  BNDRY  to  re-calculate  the  variable  boundary  data.  Note 
that  both  predicted  and  corrected  values  of  the  source  term 
vector  (S  and  SB  respectively)  are  used  in  the  corrector  step, 
which  requires  storage  of  the  S  array. 


C - CLCULATION  OF  TINESTEP,DT,  USING  THE  COUSANT-F-L  CONDITION: 

ENX=XD/RQ 

IAAX=2.»£NX 

CX=tMAX+l./EMO 

CR=1.+1./EX0 

DTXX=CFL»DX/CX 

DREF=DRR(JRM1) 

DTRR=CFL»DR£F/CR 
DT=AM1N1 (DTXX,  DTRfi) 

C - MBR  VERSION,  MACCORMACX  INTEBRflTIDN  OF  COLDFLO  .  12/5/83 

C.....FQR  SDATA1.  NONUNIFORM  DR-MESH.  NOTE:  DRR(J)=INPUT  VECTOR. 
SUH4L 

DO  1  J=l, JRMl 

1  SUM=SUN+DRR(J) 

RR(1)=6. 

DO  2  J=2, JRMl 
I=J-1 

DRR(I)=DRR(I)/SUM 

RR(J)=RR(1)+DRR(I) 

R3(J)=(DRR(I)+DRRiJ/}*RR(J)/2. 

2  R2(I)=RR(J)«2-RR(I)«2 
RR ( JRR ) =RR ( JRMl ) +DRR ( JRH 1 ) 

R2(JRR)=RR(  JRR)  **2-RR  ( JRMl )  «2 
RRS= (RR ( JRR) -RR ( JRR-2) ) /DRR (JRMl) 

RSM1=RRS-1. 

C - UNIFORM  X-MESH: 

xxu>=«. 

DO  3  K=2,KXX 

3  XX (K)=(K-1)*DX 

C - INITIAL  DATA — PRIMITIVE  VARIABLES  FIRST:DENSITY,  RAD. VEl, 

C  AXIAL  VEL,  PRESSURE.  FROM  CHICK'S  (1966!  ANALYSIS. 

PI=4.*ATANll. ) 

PIT=PI/2. 

VR(1)=0. 

VX(1!=1. 

DO  4  J=2,JRR 
AR6=PIT*RR(J)##2 
VR < J) =— SIN(ARS) /RR(J) 

4  VX(J)=COS(ARG) 

VX(JRR)=8. 

DO  S  K=1,KXX 
ARB=PI*XX(K) 

ARS2=AR6»2 

ARG3=1./(1.+(1.+AR62)«6AM3T) 

AR63=1.-6AM3T«(1.+ARG2) 

DO  5  J=l, JRR 
VIC(1,  J,K)=I. 

VIC(2,J,K)=VRU) 

VIC (3, J,  K)=VX ( J) *ARfi 

5  V1CI4,  J,K)=ARG3 
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*■ - TEWPORftRY  IfPUT:  THE  WALL  ENTHALPIES, AT  X=®,  HWO.  ft!  INI, HUP 

HUM. 

HMP=HUO 

C - AT  X=«,  M  RADIAL  VELOCITY=0: 

DO  51  J=2,JAH1 
51  VIC(2,J,1>=0. 

C - THE  INITIAL  DATA— ACTUAL  DEPENDENT  VARIABLES: 

DO  6  J-l.JRK 
DO  6  K=1,KH 
RHO=VIC(l,J,K) 

UB(1,  J,K)=RH0 
UB(2,  J,K)=RHD*VIC(2, 

UB(3,J,K)=RHD*VIC(3,J,K) 

6  UB(4,  J,K)=VIC(4,  J,K) 

DO  7  >1,  JRR 

DO  7  K=1,KXX 
DO  7  H=1,HXX 

7  U(H,J,K)=UB(N,J,K) 

C. . NOTE:  THE  FOREGOING  INCLUDES  THE  BOUNDARY  DATA.  AlL  OF  THE 

C  ABOVE  SOLUTIONS  SATISFY  THE  INITIAL  SET  OF  BCS. 

WRITE  (6, 1008)  RSTAR,  RO,  XSTAR,  XO,  VZERO,  VO,  PSTAR,  PO 
TSCAL£=RSTAR/VZERO 

WRITE (6, 1005) RHOG,  SSND, TSTAR, RUJ,  WEAR, HSTAR, CP, 

*  REO, IPRNT1, PRN, ITHAX, EHO, TSCALE, GAMA, RSTAR, CFL, DREF,  DX, DT 
1000  FORMAT (1H1,BX, 'MOTOR  RADIUS («)=’ ,F5. 2, 10X,’RO (DIMENSIONLESS)  =' 

*  ,F5. 2/9X, ’ MOTOR  L£NGTHIN)=‘,F5.2,12X,’X0(DINENSI0NLESS)=’,F5.2/ 

|  •  *  ,9X,’ INJECTION  VELOCITY <N/SEC)=’,F5. 2,’  VQ(DI*NSIQNL£SS)=’,F5.2 

*  ,/9X, 'GAS  PRESSURE (N/N«2)s' ,  1PE9. 2, 3X, ’PO(DIMENSIQNLESS)a' , 

*  0PF5.2) 

1005  FORMAT (9X, '  GAS  DENSITY <KB/0H«3 )  =» , F6. 3/ 

«  9X,  'SPEED  OF  SOUND  (M/SEC)  =\1PE9. 2/ 

*  9X,  ’  TSTAR=' , 0PF7.  3,  /, 9X, '  RUJ-' , F6. 3,  /, 9X, 

*  ’  WBAR-’ ,  F6. 3,  /,  9X, '  HSTAR=' ,  1PE9. 2,  /,  9X, 

*  'CP=',  1PE9.2,/, 

*  9X, ’RE0=’,1PE9. 2, BX,’ PRINTOUT  EACH=',  15,’  TIMES TEPS1 

*,  /9X,  ’  PRN=’ ,  0PF6. 3, 1 IX,  ’  TOTAL  RUN  DURATIONS ,  15,  ’  TIHESTEPS* 

*,  /9X,  ’  Ef10=’ ,  1PE9. 2,  BX,  ’  TSCALE  (SEC)  =' ,  1PE9. 2 
«,/9X,  ’  BAHft=’ , »PF6. 3, 10X,  ’  XSCALE (M) =’,  1PE9.2 
*,  /9X, '  CFl='  ,  0PF6. 3/9X, » DREF=’ ,  F7. 4/9X,  ’DX=’ ,  F7. 4/9X,  ’  DT=’ ,  1PE9. 

*  2 ,////) 

D1FPR=.TRUE. 

TIME*. 

L=€ 

CALL  PRINT2(TIME,L) 

DIFFPR=. FALSE. 

RETURN 

END 


following  the  printout  of  the  constant  parameters  at  the  last 
section  of  Subroutine  SDATA ,  subroutine  PRINT  is  called  to  dump 
the  initial  profiles,  in  the  form  of  X-maps  and  R-maps. 


3^6 _ PRINTOUT 

3.6.1 _ Input  Data  and  Parameters  Calculated  in  1 SDADA 


C - HDSC05  OUTPUT  FILE/3.25. 1964 — 

MOTOR  RADIUS (M)  =  8.05  RO(DIMENSIDNLESS)*  1.00 

NOTOR  LENGTH (M)=  0.55  XO(DINENSIQK£SS>=11.00 

INJECTION  VELOCITY (N/SEC)*  1.00  VO(DIMENSIONL£SS>=  1.01 

BBS  PRESSURE (N/NH2)=  2.WE+05  POIDINENSIONLESS)*  1.00 

GAS  DENSITY (RS/N»*3)=  2.423 

SPEED  OF  SOUND  (N/SEC)*  3.40E+02 

TSTAR=278.000 

(UJ=  8.314 

4IBAR*  0.028 

USTAR=  2. 89E+05 

CP*  1.04E+03 

RE&=  1.21E+04  PRINTOUT  EAC*=  200  TIMESTEPS 

PRN*  0.626  TOTAL  RUN  DURATION*  1201TINESTEPS 

EMO=  2.94E-03  TSCALE (SEC)  =  5.00E-02 

GAMA*  1.400  X  SCALE  (M)=  5.00E-02 

CFL*  0.500 

DREF*  0.0850 

DX=  1.0000 

DT=  7.33E-46 
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3.6.2  Initial  Data  (Culick  Profiles)  for  Datum  Case 
TIK  *  0.0  TIKSTEP  NO.  =  0 


X-NAP,  AT  DISTINCT  RADIAL  POSITIONS 


RHO  =  DIN. LESS  DENSITY 


Y2 

XZ=  0.0 

XZ=  0.2000 

XZ=  0.5500 

XZ=  0.8000 

XZ=  0. 9b00 

XZ=  0.9950 

0.0 

1.0B00E+00 

1.0000E+00 

t  (kWior.nrt 
At  wRwC'to 

1.0000E+00 

itRXWC  TWp 

1.0000E400 

1.000 

1.0000E+80 

1.0000E+00 

1.0000E+00 

1.0000E+00 

1.0000E+00 

1.0000E400 

2.000 

1  AAAATiiU 

A.  WWt+W 

4  aa(Mrjwi| 

itwncnw 

1.0000E400 

1.00O0E+00 

1.0000E400 

1000 

I.WWctW 

t  ACMiACxiaa 
it 

1.0000E+00 

1. 0000E+00 

1.0000E+00 

1.0000E+00 

4.000 

I  AAAATiOd 

UODOvCmw 

f  (WXW.ftA 

itOOWCnw 

i  ftViPrr ifv 
At  BWCCnw 

1.0000E+00 

At  WWtMW 

1.0000E+00 

5.000 

1  AAMkCxAA 

i  ftflaflCxiki 
It 

4  fl(vwtr  i  na 
it  vmOkTw 

1. 0000E+00 

1.0000E+00 

1.0000E+00 

6.000 

4  ftftMr.M 

KPBOTCTCHD 

1  ■  IMk 

It  PWW.’TC 

1  jflflflr  a  /wo 
It  WHPDfc+W 

1.0000E+00 

1.0000E400 

1.0000E+00 

7.000 

1.0000E+00 

1  ,  /Ml 

At  NWH.YW 

\  QOMCaI M> 

1.0000E+00 

1.0000E+00 

1.00B0E+W 

8.000 

4  OAQATilU 
it  WWW. 1  w 

«  (Mwmr .  (Vi 

it  WWt " Wu 

At  VWtl'W 

1.0000E+00 

1.0000E+00 

At  IwvL"  w 

9.000 

1.0000E+0B 

1.00006*00 

1.0000E+00 

1.0000E+08 

1.0000E+00 

1.0000E+00 

10.000 

4  (www  ■  nsk 
it  WW.'W 

4  (Wlflflr  ,  A.1 
it  *  w 

1.0000E+00 

1.0000E+00 

1.0000E+00 

1.B000E+00 

11.000 

i  AioariAA 
it  UVWCnw 

4  OAKW./WI 

At  UVWUiTlw 

4  AQ/WhCxAQ 

At  wwltW 

1.0000E+00 

1.0000E+00 

1.0000E+00 

VR  *  DIM. LESS  RADIAL  VELOCITY 

YZ 

xz=  0.0 

XZ=  0.2000 

XZ=  0.5500 

XZ=  0.B000 

XZ=  0.9b00 

XZ=  0.9950 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

1.000 

0.0 

-11395E-01  -fl.  3179E-01 

-1.0554E+00 

-1.0338E+00 

-1.0049E+00 

2.000 

0.0 

-3. 1395E-01  -B. 3179E-01 

-1.0554E+00 

-1.0338E+00 

-1.0049E+B0 

1000 

0.0 

-3. 1395E-01  -6.3179E-01 

-1.0554E+00 

-1.0338E+00 

-1.0049E+00 

4.000 

0.0 

-3. 1395E-01  -fl.  3179E-01 

-1.0554E+00 

-1.0338E+00 

-1.0049E+00 

5.000 

0.0 

-1 1395E-01  -B.3179E-01 

-1.0554E+00 

-1.0338E+00 

-1.0049E+00 

6.000 

0.0 

-11395E-01  -B.  3179E-01 

-1.0554E+00 

-1.0338E+00 

-1.0049E+00 

7.000 

0.0 

-3. 1395E-01 

-B. 3179E-01 

-1.0554E+00 

-1.033BE400 

-1.0049E+00 

B.000 

0.0 

-3. 1395E-01  -6.3179E-01 

-1.0554E+00 

-1.0336E+00 

-1. 0049E+08 

9.000 

0.0 

-3.1395E-01  -B.3179E-01 

-1.0554E+00 

-1.0338E+00 

-1.0049E+00 

10.000 

0.0 

-3. 1395E-01  -6.3179E-01 

-1.0554E+00 

-1.0338E+00 

-1.0049E+00 

11.000 

0.0 

-3. 1395E-01 

-8.3179E-01 

-1.0554E+0B 

-1.0338E+B0 

-1.0049E400 

VX  =  DIN. LESS  AXIAL  VELOCITY 

YZ 

XZ=  0.0 

XZ=  0.2000 

XI-  0.5500 

XZ=  0.8000 

XZ*  0. 9b00 

XI-  0.9950 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

1.000 

3. 1416E+00 

11354E+00 

2.7936E+00 

1.6833E+00 

3. 8591E-01 

4.9227E-02 

2.000 

6.2832E+00 

6.27BBE+00 

5.5871E+00 

3.3667E+00 

7.71B3E-01 

9.  B454E-D2 

1000 

9.4248E+00 

9. 4062E+00 

8. 3807E+00 

5. 0501E+00 

1.1577E+00 

1.4768E-01 

4.000 

1.2566E+01 

1.2542E+01 

1.1174E+01 

6.7334E+00 

1.5437E+00 

1.9691E-01 

5.000 

1.570BE+01 

1.5677E+01 

1.3968E+01 

8.4168E+08 

1.9296E+B0 

2.4614E-01 

6.000 

1.6850E+01 

1.8812E+01 

1.6761E+01 

1.0100E+01 

2.3155E+00 

2.9536E-01 

7.000 

2. 1991 E+01 

2. 1948E+01 

1.9555E+01 

1.1783E+81 

2.7814E+00 

3.4459E-01 

B.000 

2.5133E+01 

2.50B3E+01 

2.2348E+01 

1.3467E+01 

3. 0873E+00 

3. 9382E-01 

9.000 

2.B274E+01 

2.B219E+01 

2.5142E+01 

1.5150E+01 

3.4732E+00 

4. 4304E-01 

10.000 

3. 1416E+01 

11354E+01 

2. 7936E+01 

1.6833E*01 

18591E+B0 

4. 9227E-B1 

11.000 

3.4558E+01 

3.4489E+01 

3.0729E+01 

1.8517E+01 

4.2450E+00 

5.4150E-01 

Continued 
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P  *  DIM. LESS  STATIC  PRESSURE 
YZ  XZ=  0.0  XZ=  0.2000 
0.0  9.9999E-01  9. 99996-01 

I. 000  9.9993E-01  9.9993E-01 
2.000  9.9975E-01  9.9975E-01 
3.000  9.9946E-01  9.9946E-01 
4.000  9.9904E-01  9.9904E-01 
5.000  9.9B50E-01  9.9850E-01 
6.000  9.9764E-01  9.9784E-01 
7.000  9.9706E-01  9.9706E-01 
6.000  9.9617E-01  9. 9617E-01 
9.000  9.961SE-01  9.9515E-01 
10.000  9.9402E-01  9.9402E-01 

II. 000  9.9276E-01  9.9276E-01 


XZ=  0.5S00  XZ=  0.6000  XZ=  0.9608 
9.9999E-01  9. 9999E-01  9.9999E-01 
9.9993E-01  9. 9993E-01  9.9993E-01 
9.9975E-01  9. 9975E-01  9.9975E-01 
9.9946E-01  9.9946E-01  9.9946E-01 
9.9904E-01  9.9904E-01  9.9904E-01 
9. 9650E-01  9.9850E-01  9.9850E-01 
9.9784E-01  9.9764E-01  9.9764E-01 
9.9706E-01  9.9706E-01  9.9706E-01 
9.9617E-01  9.9617E-01  9.9617E-01 
9.9515E-01  9.9515E-01  9.9515E-01 
9. 9402E-01  9.9402E-01  9.9402E-01 
9.9276E-01  9.9276E-01  9.9276E-01 


XZ=  0.9958 
9.9999E-01 
9.9993E-01 
9. 9975E-01 
9.9946E-01 
9. 9904E-01 
9. 9850E-01 
9.9784E-01 
9.9706E-01 
9.9617E-01 
9.9515E-01 
9. 9402E-01 
9.9276E-01 


R-WP,  AT  DISTINCT  AXIAL  POSITIONS 


*• 


RHO  =  DIM. LESS  DENSITY 


YZ 

XZ=  0.0 

XZ=  2.0000 

XZ-  4.0000 

XZ:  6.0000 

XZ=  8.0000 

XZ=10. 0000 

0.0 

!•  wmcnv 

1  ********** 
la 

la  WWtMW 

t  flfiw  i  nn 

la 

1.0000E+00 

0.050 

i  flwiOTr  i  fin 

la  WwCmw 

la  UUUULTTO 

1.0000E400 

4  OMAriAA 

liVMtnV 

la  WW0C.TXW 

1.0000E+00 

0.200 

1.0000E+00 

I.IWwCtw 

la 

<  rmtftr ,  r%n 
la  IMRUL 1  w 

1.0000E+00 

1  4wv/v%r 
la  WWlTw 

0.400 

1.0000E+00 

4  OMATiAA 

1.WW.IW 

i  Mkkiciika 

la  WWtMW 

1.0000E+00 

1.0000E+00 

1.0000E+00 

0.550 

la  owwcnw 

|  MfW  iiM 

laWWttW 

1.0000E+00 

la 

1.0000E+00 

0.700 

1. 0000E+00 

i  anaaciAA 

4  aQQACiaa 

1.0000E+00 

1.0000E+00 

1.0000E+00 

0.800 

i  nycifla 

ItVWVC'w 

a  aaanr.M 

laUWWL^W 

1.0000E+00 

1.0000E+00 

I.0000E+00 

1.0000E+00 

0.9B0 

1.0000E+00 

la 

1.0000E+00 

1.0000E+00 

1. 0000E+00 

1.0000E+00 

0.960 

1  AkMCxlMk 
1.  UWVL 1 W 

laVWVCTw 

1.0000E+00 

1.0000E+00 

1.0000E+00 

iaWWl'lW 

0.985 

1  ftAOOCxAa 

la  «wwtTw 

la  VPflPCTVO 

1.0000E+00 

1.0000E+00 

1.0000E+00 

1.0000E+00 

0.995 

laVMRwC^W 

1. 0000E+00 

1.0000E40 

1.0000E400 

1.0000E+08 

1.000 

i  aaaocijtt 

i  ttkoocxAa 

la  wWlTW 

1.0000E+00 

1.00O0E+00 

1.0000E+00 

1.0000E+00 

VR  =  DHL  LESS  RADIAL  VELOCITY 

YZ 

XZ=  0.0 

XZ-  2.0000 

n-  4.0000 

XZ=  6.0000 

XZ=  8.0000 

XZ=10. 0000 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.050 

0.200 

0.400 

0.550 

0.700 

0.600 

0.900 

0.960 

0.9B5 

0.995 

1.000 


0.0 

-7.B539E-02 

0.0 

-3. 1395E-01 

0.0 

-6.2172E-01 

0.0 

-8.3179E-01 

0.0 

-9.9416E-01 

0.0 

-1.0554E400 

0.0 

-1.0620E+00 

0.0 

-1.0338E+00 

0.0 

-1.0141E+00 

0.0 

-1.0049E+00 

“1.0000E+00 

-1.0000E+00 

-7.8539E-02  -7.B539E-02 
-3. 1395E-01  -3. 1395E-01 
-6.2172E-01  -6.2172E-0I 
-B.3179E-01  -6.3179E-01 
-9.94I6E-01  -9.9416E-01 
-1.0554E+00  -1.0554E+00 
-1.0620E+00  -1.0620E+00 
-1.0338E+00  -1.0338E+00 
-1.0141E+00  -1.0141E+00 
-1.0049E+00  -1.0049E+00 

IMAATiOa  4  AMATiM 
a  1 W  1  a  Wwfc  *  V%l 


-7.6539E-02  -7.6539E-02 
-3. 1395E-01  -3. 1395E-01 
-6.2172E-01  -6.2172E-01 
-B.3179E-01  -8. 3179E-01 
-9.94I6E-01  -9.9416E-01 
-1.0554E+00  -1.0554E+00 
-I.0620E+00  -1.0620E+08 
-1.0338E+00  -1.0338E+00 
-1.014IE+00  -1.0141E+00 
-1.0049E+00  -1.0049E+00 
-1. 0000E+W  -1.0000E40 


Continued 
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VX  =  DIN. LESS  AXIAL  VELDC1TY 


Y2 

XZ=  8.0 

XZ=  2.8888 

XZ:  4.8888 

XZ:  6.0000 

XZ:  8. 0888 

XZ=10.0888 

it 

0.0 

S.2832E+80 

1.2566E+01 

1.BB50E+01 

2.5133E+01 

3. 1418E+81 

tt.060 

0.0 

6.2831E+08 

1.2566E+81 

1.8849E+01 

2.5133E+01 

3. 1416E+81 

1.280 

0.0 

6.2708E+88 

1.2542E+01 

1.8B12E+01 

2.5883E+01 

3. 1354E+81 

0.400 

0.0 

8.0858E+80 

1.2172E+81 

1.8257E+01 

2. 4343E+81 

3.8429E+81 

1.558 

0.0 

5.5871E+88 

1.1174E+01 

1.6761E+01 

2.2348E+81 

2. 7936E+01 

8.708 

0.8 

4.5121E+00 

9.0242E408 

1.3536E+81 

1.8848E+81 

2.2561E+01 

8.888 

8.0 

3.3667E+08 

6.7334E+08 

1.0100E+01 

1.3467E+01 

1.6833E+01 

0.988 

0.0 

1.B475E+08 

3.6958E+08 

5.5426E+08 

7.3901E+88 

9. 2376E+88 

0.988 

0.0 

7.71B3E-01 

1.5437E+00 

2.3155E+00 

3.0673E+88 

3.8591E+08 

0.985 

0.0 

2. 9377E-81 

5.8753E-01 

8.8138E-01 

1.1751E+00 

1.4688E+88 

0.995 

0.0 

9.B454E-02 

1.9691E-81 

2.9536E-01 

3.9382E-01 

4.9227E-81 

1.008 

0.0 

0.0 

0.0 

0.0 

8.8 

0.0 

P  »  DIN.  LESS  STATIC  PRESSURE 

YZ 

XZ =  0.8 

XZ*  2.8888 

XZ:  4.8008 

XZ:  6.0000 

XZ:  8.0808 

XZ:18.0888 

0.0 

9.9999E-81 

9.9975E-01 

9.9984E-81 

9.9784E-01 

9.9617E-01 

9.9402E-81 

0.058 

9.9999E-01 

9. 9975E-81 

9. 9904E-81 

9. 9784E-01 

9.9617E-01 

9. 9402E-01 

0.208 

9.9999E-81 

9.9975E-81 

9.9904E-01 

9.97B4E-01 

9.9617E-01 

9. 9482E-81 

0.480 

9.9999E-01 

9.9975E-01 

9.9904E-01 

9.9784E-01 

9.9617E-81 

9.9402E-81 

0.558 

9.9999E-81 

9.9975E-01 

9.9904E-01 

9. 9784E-81 

9.9617E-01 

9.9482E-81 

0.788 

9.9999E-01 

9.9975E-01 

9. 9904E-01 

9.9784E-01 

9.9617E-81 

9.9482E-01 

8.888 

9.9999E-01 

9.9975E-81 

9.9904E-01 

9.97B4E-01 

9.9617E-01 

9.9482E-01 

0.908 

9.9999E-01 

9.9975E-81 

9.9984E-81 

9.9784E-81 

9.9bl7£-01 

9.9482E-81 

8.968 

9.9999E-01 

9.9975E-01 

9.9904E-81 

9. 9784E-01 

9.9617E-81 

9.9482E-81 

0.985 

9.9999E-81 

9.9975E-81 

9.9904E-01 

9. 97846*01 

9.9617E-81 

9. 9482E-81 

0.995 

9.9999E-81 

9.9975E-01 

9.9984E-01 

9.9784E-81 

9.%17E-81 

9.9402E-01 

1.888 

9. 9999E-01 

9. 9975E-01 

9.9984E-81 

9.9784E-81 

9.9617E-01 

9. 9482E-81 

INTERMEDIATE  STEPS  PRINTOUT  FOLLOWS 


(FOR  n  <  ITMAX) 
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3^6 l2 _ Final _ Step  Printout,  Datum  Case _ (Step  No.  1201) 


T1K  *  A.  7S61E-43  TIlBTt.  .....  -1201 


HAP,  AT  DISTINCT  RADIAL  POSITIONS 


RHO  *  DIM. LESS  DENSITY 


Y2 

X2=  0.0 

XZ=  0.2000  XI-  0.5500 

0.0 

9. 8137E-01 

9.B133E-01  9.8133E-01 

1.000 

9. 8850E-01 

9. 8851E-01  9.8855E-01 

2.000 

9.8853E-01 

9.6853E-01  9.8887E-01 

3.000 

9.8653E-01 

9.8852E-01  9.8857E-01 

4.000 

9.8856E-01 

9.8852E-01  9.8885E-01 

5.000 

9. 8858E-01 

9.8b53E-01  9.8862E-0; 

s.000 

9.BE58E-01 

9.8653E-01  9.8860E-01 

7.000 

9.8858E-01 

9.  B853E-01  9.8662E-0; 

8.000 

9.8658E-01 

9.8854E-01  9.8b6JE-01 

9.000 

9.B860E-01 

9. 8659E-01  9.8888E-01 

10.000 

9.8685E-01 

9. B889E-01  9.8676E-01 

11.000 

9.8871E-01 

9.B8B0E-01  9.B889E-01 

VR  =  DIN.  LESS  RADIAL  VELOCITY 

YZ 

X2=  0.0 

X2=  0.2000  XZ=  0.550* 

0.0 

0.0 

0.0  0.0 

1.000 

0.0 

-2. 4673E-01  -9.7828E-01 

2.000 

0.0 

-2.4820E-01  -9.7119E-01 

3.000 

0.0 

-2.4338E-01  -9.7252E-01 

4.000 

0.0 

-2.4487E-01  -9.6701E-01 

5.000 

0.0 

-2. 4837E-01  -9.7014E-01 

8.000 

0.0 

-2. 5339E-01  -9.7080E-01 

7.000 

0.0 

-2.4248E-01  -9.7718E-01 

8.000 

0.0 

-2.5518E-01  -9.8039E-01 

9.000 

0.0 

-2.4381E-01  -9. 7789E-01 

10.000 

0.0 

-2.5297E-01  -9. 6403E-01 

11.000 

0.0 

-2.6213E-01  -9.5017E-01 

Xl=  0.8*00 
9. 8135E-01 
9.8859E-01 
9.886tE-01 
9.8c82E-01 
9.B6tbE-01 
9.8b83E-01 
9.8883E-01 
9. 8b7«E-e: 
9.867*€-0: 
9.8679E-61 
9.8688E-01 
9.8t>93E-01 


X2=  0.8*00 

0.0 

-1. 1914E+W 
-1. 1660E+08 
-1.1697E+W 
-1.  1661E+00 
-1.1756E+W 
-1. 1627E+00 
-1.1691E+00 
-1. 1853E+00 
-1. 16S7E+00 
-1. 1675E+00 
-1.1653E+00 


XI-  0. 9b00 
9. SlBlE-01 
9.86&4E-0: 
9.8889E-01 
9. 8o87E-01 
9.8878E-01 
9. 8e>9*E-0: 
9.6899E-01 
9.8t,73E-01 
9. 8701E-01 
9.87I7E-01 
9.8717E-01 
9.8/17E-01 


XZ=  0.9*5* 
9.8196E-01 
9.8^3E-0i 
9.8461E-01 
9.6".3/E-0l 
9. B369E-01 
9.B373E-01 
9. B338E-01 
9. 6c:3t,E-0i 
9.8218E-01 
9.8185E-01 
9.6077E-01 
9.7989E-01 


XZ=  0.9b0*  XZ=  0.9950 
0.0  0.0 
-9.0368E-0:  -1.003bE+00 
-1.0278E+00  -1.0188E+00 
-9.6994E-0:  -1.0l3bE+00 
-1.0551E+00  -1.0237E+00 
-1.0231E+W  -1.B199E+00 
-1.0118E+00  -1.8192E+00 
-1.0146E+0*  -1.0206E+00 
-1.0057E+0B  -1.0183E+00 
-9. 4666E-01  -1.0129E+00 
-1.0771E+00  -1.0305E+00 
-1.2078E+00  -1.048IE+W 


VX  =  DIM. LESS  AXIAL  VELOCITY 


Y2 

11=  0.0 

11=  0.2000 

XZ=  0.5500 

X2=  0.8000 

X2S  0.9b00 

XZ=  0.9950 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

1.000 

3. 1669E+00 

3. 1613E+00 

2.8102E+00 

1.6932E+00 

3.8/16E-01 

3.9751E-02 

2.000 

6.2929E+00 

8. 2794E+00 

5. 5649E+00 

3.3570E+00 

7.5780E-01 

8.6714E-02 

3.000 

9. 4458E+00 

9.4238E+00 

8. 3856E+00 

5.0507E+00 

1.1542E+00 

1.4o9*£-01 

4.000 

1.2574E+01 

1.2547E+01 

1.1161E+01 

6.7178E+00 

1.5228E+00 

1.8374E-01 

5.000 

1.5716E+01 

1.5682E+01 

1.3951E+01 

8. 3326E+00 

1.9K80E+00 

2. 351 0£ -01 

8.000 

1.8845E+01 

1.8805E+01 

1.8732E+01 

1.0075E+01 

2. 2930E+M 

2. 8389E-01 

7.000 

2. 1987E+01 

2. 1938E+01 

1.9519E+01 

1. 17522+01 

2.6719E+00 

3. 2tt4AE-01 

B.000 

2.5115E+01 

2.5057E+01 

2.2291E+01 

1.3417E+01 

3. 0407E+00 

3.6721E-01 

9.000 

2.8243E+01 

2.8180E+01 

2.5073E+01 

1.5088E+01 

3.410*2+0* 

4.0293E-01 

10.000 

3. 1348E+01 

3. 1274E+01 

2.7820E+01 

1.6717E+01 

3.7468E+00 

4.2128E-01 

11.000 

3.4452E+01 

3.4387E+01 

3. 0588E+01 

1.B345E+01 

4.0635E+0* 

4.3959E-01 
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P  =  DIPL  LESS  STATIC  PRESSURE 
¥Z  XZ=  8.0  XZ=0.20W 

0.0  9.8137E-01  9.8133E-01 

1.000  9.8125E-01  9.8123E-01 
2.000  9.8114E-01  9.8113E-01 
3.000  9.8885E-01  9.8883E-01 
4.000  9.8047E-01  9.8044E-01 
5.000  9.7998E-01  9.7995E-01 
6.000  9. 7935E— 01  9.7934E-01 
7.000  9. 7866E-01  9.7861E-01 
8.000  9.77B2E-01  9.7780E-01 
9.000  9. 7693E-01  9.7693E-01 
10.000  9.7596E-01  9.7596E-01 
11.008  9.7499E-01  9.7508E-01 


KZ=  0.5500  XZ=  0.8000  XZ=  0.9c>00 
9.8133E-01  9.8135E-01  9.8181E-01 
9.8128£-01  9.8132E-01  9.6172E-01 
9.8124E-01  9.8129E-01  9. 8164EH01 
9.8088E-01  9.8092E-01  9.8133E-01 
9.8052E-01  9. 8059E-01  9.8076E-01 
9.8802E-01  9. 8805E-01  9.8847E-01 
9. 7938E-01  9.7944E-01  9.7996E-01 
9. 7870E-01  9. 7675E-01  9.7686E-0] 
9. 7787E-0 1  9.7791E-01  9.7641E-01 
9. 7702E-01  9. 7704E-01  9.7767E-01 
9. 7605E-01  9.7607E-01  9.7659E-01 
9. 7507E-01  9.751*1-01  9.7558E-01 


XZ=  8.9950 
9.8198E-01 
9.8193E-01 
9. 8189E-01 
9.6l5'9E-01 
9.6094E-01 
9. 8078E-01 
9. 8829E-01 
9.7B96E-01 
9. 7877E-01 
9.7alB£-01 
9. 7693E-01 
9. 7575E-01 


DIFFERENCES  BETWEEN  CURRENT  TIftXrEP  NO.  1201 
AND  THE  INITIAL  DATA  U<«, J,K)-UINITIAL(«,  J,K>  : 


HUP,  AT  DISTINCT  RADIAL  POSITIONS 


RHO  =  DIN. LESS  DENSITY 

YZ  XZ=0.0  XZ=  0.2000  XZ=  0.5500  XZ=  0.8U08  XZ=  0.9bW  XZ=  0.9950 
0.0  -1.8634E-02  -1.8672E-02  -1.8675E-02  -1.8647E-02  -1.8193E-02  -1.W21E-02 
1.000  -1.3495E-02  -1.3494E-02  -1.3447E-02  -1.3410E-02  -1.3165E-0E  -1.5371E-02 
2.000  -1.3474E-02  -1.3468E-02  -1.3333E-82  -1.3337E-02  -1.310faE-02  -1.5383E-02 
3.000  -1.3466E-02  -1.3483E-02  -1.3438E-02  -1.3379E-02  -1.3129E-02  -1.5633E -0c' 
4.000  -1. 3424E-02  -1.3482E-02  -1.3354E-02  -1.3339E-02  -1.3239E-02  -1.6114E-6£ 
5.000  -1.3416E-02  -1.3472E-02  -1.3377E-B2  -1.3368E-02  -1.3103E-02  -1.626bE-te 
6.000  -1.3436E-02  -1.3469E-02  -1.34&4E-02  -1.3372E-02  -1.3008E-02  -1.6644E-02 
7.000  -1.3423E-02  -1.3471E-02  -1.3381E-02  -1.3304E-02  -1.3E73E-02  -1.7637E-02 
B.000  -1.3425E-02  -1.3462E-02  -1.3368E-82  -1.3303E-02  -1.2994E-02  -1.7821E-02 
9.000  -1. 3402E-02  -1.3412E-02  -1.3321E-02  -1.3212E-fc  -1.2627E-02  -1.8349E-02 
10.000  -1.3346E-02  -1.3306E-02  -1.3215E-02  -1.3143E-02  -1.2828E-02  -1.9231E-02 
11.000  -1.3291E-02  -1.3201E-02  -1.3H0E-82  -1.3873E-02  -1.2838E-02  -2.0114E-02 


VR  =  DIN.  LESS  RADIAL  VELOCITY 

XZ=  0.2000  XZ=  0.5500  XZ=  0.8800  XZ=  0.9b08  XZ=  0.9950 

0.0  0.0  0.0  0.0  0.0 

6.721BE-02  -1.4649E-01  -1.3595E-01  1.2990E-01  1.27B9E-03 
6.5748E-02  -1.3940E-01  -1.1056E-01  6.0110E-03  -1.3935E-82 
7.0591E-02  -1.4072E-01  -1.1428E-01  6.38342-02  -8.6632E-03 
6.90B3E-82  -1.3522E-01  -1.1065E-01  -2.1356E-02  -1.B764E-02 
6. 5582E-02  -1.3635E-01  -1.2044E-01  1.0669E-02  -1.5805E-08 
6.0562E-02  -1.3901E-01  -1.0731E-01  2.1954E-02  -1.4344E-02 
7. 1472E-02  -1.4539E-01  -1.1369E-01  1.9l44E-fc  -1.5697E-02 
5. 8774E-02  -1.4660E-01  -1.0994E-01  2.8833E-02  -1.3985E-02 
7.0144E-02  -1.4610E-01  -1.1433E-01  8.7H6E-02  -7.9870E-03 
6.0981E-02  -1.3224E-01  -1.1214E-01  -4.3362E-02  -2.5578E-B2 
5.1820E-02  -1. 1B38E-01  -1.0994E-01  -1.7384E-01  -4.320^-08 


YZ 

XZ=  0.0 

0.0 

0.0 

1.000 

0.0 

2.000 

0.0 

3.000 

0.0 

4.000 

0.0 

5.000 

0.0 

6.000 

0.0 

7.000 

0.0 

8.000 

0.0 

9.000 

0.0 

10.000 

0.0 

11.000 

0.0 
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VX  «  DIB. LESS  AXIAL  VELOCITY 

YZ  XZ=  0.0  XZ=  0.2000  XZ=  0.5500  XZ=  0.0000  XZ=  0.9b00  XZ=  0.9950 

0.0  0.0  0.0  0.0  0.0  0.0  0.0 

I. 000  8.04696-03  8.23816-03  5.29696-03  3.14286-03  4.0304E-04  -3.0165E-03 

2.000  1.54306-03  1.3724E-03  -3. 58366-04  -1.5509E-03  -2. 23206-03  -1.5502E-03 
3.000  2.23206-03  1.06766-03  5.21126-04  7.23496-05  -3.71066-04  -8.2405E-05 

4.000  6.4227E-04  3.9673E-04  -1.0210E-03  -1.23966-03  -1.6592E-03  -1.0403E-03 
5.000  5. 28266-04  2.9045E-04  -1.0530E-03  -1.52706-03  -1.3711E-03  -7.0257E-04 
6.000  -2. 33146-04  -3.70756-04  -1.5332E-03  -1.3195E-03  -1.19106-03  -6.0007E-04 
7.000  -2.09556-04  -4.38896-04  -1.6521E-03  -1.41086-03  -1.3402E-03  -7.3455E-04 
6.000  -7. 1641E-04  -1.03766-03  -2.29496-03  -1.96936-03  -1.B53BE-03  -1.65666-63 
9.000  -1.09286-03  -1.3567E-03  -2.4447E-03  -2.21336-03  -2.2342E-03  -1.41656-03 
10.000  -2. 16436-03  -2.54466-03  -3.68756-03  -3.72156 -03  -3.57686-03  -2.2598E-03 

II. 100  -3.05116-03  -3. 53596-03  -4.72106-03  -4.96266-03  -4.67396-03  -2.94896-03 


P  «  DIM. LESS  STATIC  PRESSURE 
YZ  XZ=  0.0  XZ=  0.2000 
0.0  0.0  4.22516-03 

I. 000  1.25026-02  1.51036-02 
2.000  2.50216-02  2.60386-02 
3.000  5.79896-02  6.0486E-02 
4.000  1.00466-01  1.0394E-01 
5.000  1.5614E-01  1.58826-01 
6.000  2. 26926-01  2.27326-01 
7.000  3.04266-01  3.0967E-01 
8.000  3.9891E-01  4.01686-01 
9.000  4.9924E-01  4.99716-01 
10.000  6.0954E-01  6.09136-01 

II. 000  7. 2007E-01  7. 1884E-01 


XZ=  0.5500  XZ=  0.0800 
4.59956-03  1.43736-03 
9.1191E-03  4.9404E-03 
1.36796-02  8.49736-03 
5.4590E-02  4. 99196-02 
9.4486E-02  8.68206-02 
1.5077E-01  1.48256-01 
2.23486-01  2. 16556-01 
3.00316-01  2. 9458E-01 
3.94036-01  3.89556-01 
4.89816-01  4. 87396-01 
5. 99946-01  5. 97236-01 
7. 10366-01  7.0737E-01 


XZ=  0.9bW  XZ=  0.9958 
-4. 94076-02  -6. 86726-02 
-3. 98956-02  -6.3b84E-02 
-3.03686-02  -5.86566-02 
3.86416-03  -2.5471E-02 
6.61076-02  4.63346-82 
1.01056-01  6.8129E-02 
1.56086-01  1.20886-01 
2.8178E-01  2.70676-01 
3.32676-01  2. 92156-01 
4.15756-01  3. 67526-01 
5. 38406-01  5.0007E-01 
6. 61396-01  6. 3300E-0; 


SKIN  FRICTION(CF)  t  STANTON  NO. (ST) 

X  CF  ST  U1(U)/U2(CL) 

1.00  1.60716-05  -5.25886-04  -3.20096-01 
2.00  9.08346-06  -2.6062E-04  -1.6108E-0; 
3.00  6.67586-06  -1.67656-04  -1.87316-01 
4.08  4.71136-86  -1. 24816-04  -8.06686-02 
5.0O  3. 85916-06  -9.30456-05  -6.44946-82 
6.00  3. 24106-06  -7. 32406-85  -5. 37B7E-02 
7.00  2.7547E-I6  -6.3609E-05 -4.6101E-02 
8.00  2.3604E-06  -5.0464E-05  -4.0359E-02 
9. N  2. 04826-06  -4. 23446-05  -3. 5887E-02 
11.10  1.73806-06  -3. 73326-05  -3. 2332E-82 
11.00  1.50146-06  -3. 33576-05  -2.94176-02 
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R-KAP,  AT  DISTINCT  AXIAL  POSITIONS 


RHO  *  DIN.  LESS  DENSITY 


YZ 

XZ=  0.0 

n-  2.0000 

XZ=  4.0000 

0.0 

9.8137E-01 

9.8653E-01 

9.B658E-02 

01060 

9.8137E-01 

9.8653E-01 

9.8658E-01 

0.200 

9.8133E-01 

9.8653E-01 

9.8652E-01 

0.400 

9.8111E-01 

9.8641E-01 

9.8640E-01 

0.5S0 

9.8133E-01 

9.8667E-01 

9.0665E-01 

0.700 

9.B118E-01 

9.6664E-01 

9.B655E-01 

0.800 

9.8135E-01 

9.8666E-01 

9.B666E-01 

0.900 

9.8131E-01 

9.8653E-01 

9.8657E-01 

0.960 

9.B1B1E-01 

9.B689E-01 

9.B676E-01 

0.965 

9.  B196E-01 

9.8657E-01 

9.8623E-01 

0.995 

9.819BE-01 

9.  8461E-01 

9. B3B9E-01 

1.000 

1.0088E+00 

9.8190E-01 

9.8094E-01 

VR  =  DIN. LESS  RADIAL  VELOCITY 

YZ 

XZ=  0.0 

XZ=  2.0000 

XZ=  4.0000 

0.0 

0.0 

0.0 

0.0 

0.050 

0.0 

-9. 4916E-02 

-9.4648E-02 

0.200 

0.0 

-2.4820E-01 

-2. 44B7E-01 

01400 

0.0 

-B.0947E-01 

-B.B589E-01 

0.550 

0.0 

-9.7119E-01 

-9. 6701E-01 

0.700 

0.0 

-1. 19526+08 

-1. 1947E+00 

0.600 

0.0 

-1. 1660E+00 

-1. 1661E+00 

0.900 

0.0 

-1.1461E+00 

-1.1614E+88 

0.960 

0.0 

-1.0278E+00 

-1.0551E+00 

0.965 

0.0 

-1.0227E+00 

-1.0343E+06 

0.995 

0.0 

-1.0188E+00 

-1.0237E+00 

1.000 

-1.0000E+00 

-1.0184E+00 

-1.0194E+06 

VX  =  DIN. LESS  AXIAL  VELOCITY 


YZ 

XZ=  0.0 

U-  2.00W 

XZ=  4.0000 

0.0 

0.0 

6.2929E+00 

1.2574E+02 

0.050 

0.0 

6. 2929E+80 

1.2574E+01 

0.200 

0.0 

6.2794E+00 

1.2547E+81 

0.400 

0.0 

6. 0696E+00 

1.2169E+01 

0.550 

0.0 

5.5B49E+08 

1.1161E+01 

0.700 

0.0 

4.5871E+08 

9.0119E+08 

0.B00 

0.0 

3.3570E+08 

6.7178E+00 

0.900 

0.0 

1.B37BE+08 

3. 68106+00 

0.960 

0.0 

7.5786E-01 

1.52266+00 

0.965 

0.0 

2. 8026E-01 

5.6903E-01 

0.995 

0.0 

8. 87146-02 

1.8374E-01 

1.000 

0.0 

0.0 

0.0 

XZ=  S  Ml 
9.6656E-01 
9.8656E-01 
9. 8653E-01 
9.8641E-01 
9.8660E-01 
9.8660E-01 
9.B663E-01 
9.8655E-01 
9. 8699E-01 
9.8631E-61 
9.8336E-01 
9.8038E-0I 


xz=  a.M0 

9. 8656E-01 
9.&58E-01 
9. 8654E-01 
9.8636E-01 
9.B663E-01 
9.8663E-01 
9.8670E-01 
9.8b5ZE-01 
9.6701E-01 
9.  BbWE-01 
9.B21BE-01 
9. 7678E-81 


XZ=18.WW 
9.B665E-01 
9.6b65E-0i 
9. 8669E-01 
9.8660E-01 
9. B678E-01 
9.8679E-01 
9. 8686E-01 
9.8o&iE-01 
9.B717E-01 
9.B569E-01 
9. 8077E-01 
9.7694E-B1 


XZ=10>  00W 

0.0 

-9.6113E-02 
~2. 5297E-01 
-6.0785E-01 
-9. 6403E-01 
-1.  i%9£+00 
-1.1675E+00 
-1. 1752E+00 
-1.0771E+00 
-1.045K+00 
-1.0305E+W 
-1.023bE+00 


XZ=  6.0080  XZ=  6.0008 
0.0  0.0 
-9.709bE-«2  -9.6l49£-8<: 
-2.5339E-01  -2.5516E-01 
-8.0669E-01  -fl. 1320E-01 
-9. 7080E-01  -9.8039E-01 
-1. 1980E+00  -1.2012E+00 
-1. 1627E+08  -1.1653E+00 
-1. 1378E+00  -1. 1373E+00 
-1.0118E+00  -1.0657E+80 
-I.0176E+00  -1.6155E+88 
-1.0192E+00  -1.01B9E+W 
-1.0201E+00  -1.0217E+80 


XZ=  6.WW 
1.8845E+01 
1.B845E+01 
1.B685E+01 
1.8248E+01 
1.6732E+01 
1. 35196+01 
1.007SE+01 
5.5279E+08 
2.2930E+00 
8.6322E-01 
2.8389E-01 
0.0 


XZ=  B.000O 
2.5llbE+01 
2.5115E+01 
2.5857E+01 
2.4383E+01 
2.2291E+01 
1.7997E+01 
1.3417E+02 
7. 352SE+00 
3.0407E+W 
1. 1352E+00 
3.6721E-01 
6.0 


XZ=10.0080 
3. 1348E+01 
3. 134BE+01 
3. 2274E+01 
3.0332E+01 
2. 7820E+01 
2.2*536+0! 
1.6717E+01 
9.  IM9E+00 
3.7468E+88 
1.3b41E+00 
4.212BE-01 
0.0 
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P  «  DIN.  LESS  STATIC  PRESSURE 
Y2  12=  A I  12=  2.  MM 
1.1  9.B137E-01  9.81146-01 

AC50  9.81376-01  9.81146-01 
0.2M  9.81336-01  9.81136-01 
A4M  9.81116-01  9.8894E-01 
0.550  9.81336-01  9.81246-01 
ATM  9.81186-01  9.81236-01 
ASM  9.81356-01  9.81296-01 
ASM  9.B131E-01  9.B109E-01 
A  960  9.B181E-01  9.81646-81 
A  985  9. 81986*01  9.B187E-01 
A 996  9. 8198E-01  9.81896-01 
1.M0  9.9999E-01  9.61986-01 


22=  4.08M 
9. 8047E-01 
9.8847E-01 
9.8044E-01 
9. 60266*01 
9.88526-01 
9.8051E-01 
9.80596-01 
9.80446-01 
9.80786-01 
9.80936-01 
9.80946-01 
9.80946-01 


12=  6.0000 

9.79356-01 

9.79356-01 
9.79346-01 
9.79136-01 
9. 79366-01 

9.79356-01 
9.79446-01 
9.79266-01 
9.79966-01 
9.80276-01 
9.80296-01 
9.80306-01 


22=  8.0000 
9.77826-01 
9.7782E-01 
9.77806-01 
9.77556-01 
9.77876-01 
9.77856-01 
9.77916-01 
9.77686-01 
9.7841E-01 
9.7o756-01 
9.7877E-01 
9.76786-01 


22=10.0000 
9.75966-01 
9.759bE-0! 
9.75966-01 
9.75616-01 
9.76056-01 
9. 75996-01 
9. 76076-01 
9.7591E-01 
9. 7659E-01 
9. 7b9l6-01 
9.76936-01 
9.7b94E-01 


DIFFERENCES  BETWEEN  CURRENT  TlNESTEP  NC.1201 
AND  THE  INITIAL  DATA  U(R,J,KKJINITIAL(N,J,K)  : 


R-NAP,  AT  DISTINCT  AXIAl  POSITIONS 


RriD  =  DIN.  LESS  DENSITY 

Y2  22=  A0  22=  2.80M  22=  4.0000  XZ=  6.08M  22=  8.00M  22=10.0800 
AB  -1.8634E-02  -1.34746-02  -1.34246-02  -1.34366-02  -1.34256-02  -1.33466-02 
A  050  -1.86346-02  -1.34746-02  -1.34246-02  -1.34366-02  -1.34256-02  -1.33466-02 
0.200  -1.86726-02  -1.34686-02  -1.34826-02  -1.34696-02  -1.34626-02  -1.33866-02 
A  400  -1.8887E-02  -1.3586E-02  -1.35986-02  -1.3591E-02  -1.36246-02  -1.34886-02 
0.550  -1.86756-02  -1.33336-02  -1.33546-02  -1.34046-02  -1.33666-02  -1.3215E-02 
A700  -1.88166-02  -1.3364E-02  -1.3451E-02  -1.34006-02  -1.3367E-62  -1.3212E-02 
0.800  -1.86476-02  -1.3337E-02  -1.3339E-02  -1.33726-02  -1.33036-02  -1.31436-02 
A 900  -1.86B9E-02  -1.3471E-02  -1.34386-02  -1.34526-02  -1.34306-02  -1.31996-02 
0.960  -l.B'.93E-02  -1. 31066-02  -1.32396-02  -1.30066-02  -1.29946-02  -1.28286-02 
A 985  -1.80236-02  -1.34286-02  -1.37706-02  -1.36896-02  -1.40026-02  -1.4312E-02 
0.995  -1.80216-02  -1.53896-02  -1.61146-02  -1.66446-02  -1.78216-02  -1.9231E-02 
1.000  0.0  -1.B1036-02  -1.90626-02  -1.9702E-02  -2.12226-02  -2.3064E-02 


VR  =  DIN. LESS 

yz  xz=e. 


8.e 

0.05* 
A  200 
A  480 
0.550 
A  7M 
0.800 
A908 
A960 
A965 
A995 
1.000 


0.0 

0.0 

A0 

0.0 

A0 

1.0 

1.0 

0.0 

1.0 

0.0 

0.0 

0.0 


RADIAL  VELOCITY 
0  XZ=  2.08M 
0.0 

-1.63776-02 
6.57486-02 
-1.87756-01 
-1.39486-01 
-2.01046-01 
-1. 10566-01 
-8.41566-02 
6.01186-03 
-8.57266-03 
-1.39356-02 
-1.84366-02 


22=  4.0808 

A0 

-1.61096-02 
6.98836-02 
-1.8416E-01 
-1.35226-01 
-2.08516-01 
-1. 18656-01 
-9. 93616-02 
-2.13566-02 
-2.03536-02 
-1.87646-02 
•1. 9432-02 


XZ=  6.08M 
A0 

-1.85576-02 
6.0562E-02 
-1.8497E-01 
-1.39016-01 
-2.03826-01 
-1.0731E-01 
-7.58056-02 
2. 1954E-02 
-3.43896-03 
-1.43446-02 
-2.08986-02 


22=  8.0888 

A0 

-1. 7609E-02 
5.87746-02 
-1.91486-01 
-1.48686-01 
-2.07026-01 
-1.0994E-01 
-7.5351E-02 
2.80336-02 
-1.41056-03 
-1.3S85E-02 
-2.16826-02 


22=10.0808 

A0 

-1.75736-02 

6.0981E-02 

-1.8612-01 

-1.32246-01 

-2.0271E-01 

-1.1214E-01 

-1.1319E-01 

-4.3362-02 

-3.0882-02 

-2.55786-02 

-2.36086-02 


Continued 
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VX  =  DIM. LESS  AXIAL  VELOCITY 


YZ 

XZ=  0.0 

XZ=  2.0000  XZ=  4.0008  XZ=  6.0080 

0.0 

0.0 

1.54386-03  6.42276-84  -2.3314E-04 

0.050 

0.0 

1.55086-83  6.5801E-84  -2.2504E-04 

0.200 

0.0 

1. 3724E-03  3.98736-04  -3.70756-04 

0.400 

0.0 

6.05006-04  -2.09546-04  -9.33366-04 

0.550 

0.0 

-158366-04  -1.02186-03  -1.53326-03 

0.718 

0.0 

-8.8596E-04  -9.82266-04  -9.35946-84 

0.800 

0.0 

-1. 55096-03  -1.23966-03  -1.3195E-83 

0.900 

0.0 

-1.55036-03  -1.11636-03  -7.75866-04 

0.968 

0.0 

-2.23206-03  -1.65926-03  -1.19106-03 

0.965 

0.0 

-2.14926-03  -1.4725E-03  -9.5897E-04 

0.995 

0.0 

-1.5502E-03  -1.04836-03  -6.0887E-04 

xz=  e.0^  xz=ie.«t^ 

-7.1641E-04  -2. 1643E-03 
-7.88526-04  -2. 15656-03 
-1.0376E-03  -2.54466-83 
-1.58526-03  -3.88526-83 
-2.29496-03  -3.6875 6-03 
-2.03756-83  -3.41406-83 
-1.9693E-03  -3.72156-03 
-1.49156-03  -3. 04646-83 
-1. 85386-03  -3. 5760E-03 
-1.58666-03  -3.33466-03 
-1.05866-03  -2.25986-03 


1.000  0.0  0.0  0.0 


0.0  0.0  0.0 


P  *  DIM.  LESS  STATIC  PRESSURE 

YZ  XZ=  0.0  XZ=  2.0000  XZ=  4.0000  XZ=  6.0000  XZ=  8.0000  XZ=10.W08 
0.0  -1. 862BE-02  -1.B612E-02  -1.85666-02  -1.8495E-02  -1. 83466-02  -1.8054E-02 

0.050  -1.8628E-02  -1.8612E-02  -1.85666-02  -1.8495E-02  -1.83466-02  -1.80546-02 
0.200  -1.B666E-02  -1.8621E-02  -1.8597E-02  -1.84986-02  -1.8370E-02  -1.B051E-02 
0.400  -1.8881E-02  -1.88196-02  -1.87906-02  -1.87096-02  -1.8617E-02  -1. 82046-02 
0.550  -1.86696-02  -1.B511E-02  -1.85136-02  -1.8464E-02  -1.83036-02  -1.79706-02 
0.700  -1.88106-02  -1.85226-02  -1.85326-02  -1.B496E-02  -1.83176-02  -1.80246-02 
0.860  -1.86416-02  -1.84656-02  -1.84456-02  -1.B403E-02  -1.8263E-02  -1.79466-02 
0.900  -1.86836-02  -1.B663E-02  -1.85986-02  -1.85816-02  -1.649lE-08  -1.61086-03 
0.96B  -1.81876-02  -1.81166-02  -1.82686-02  -1.78836-02  -1.77596-02  -1.7427E-82 
0.965  -1.80176-02  -1.78826-02  -1.81086-82  -1.75706-82  -1.74216-82  -1.7l0bE-02 
0.995  -1.80156-02  -1.78656-02  -1.8102E-82  -1.7552E-02  -1.7400E-02  -1.70886-02 
i  1.000  0.0  -1.7658E-02  -1.B1WE-02  -1.75446-02  -1.73906-02  -1.70806-82 
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Recently  Flandro  [15]  has  carried  out  a  theoretical  analysis  for  a  burning 
pellant  in  a  cylindrical  grain,  under  the  effect  of  incident  acoustic  waves, 
etailed  formulation  was  derived  with  a  double  expansion,  in  terms  of  both 
erse  Reynolds  number  as  well  as  Mach  number  (independent  small  parameters), 
onsteady  premixed  combustion  zone  was  considered  near  the  propellant  surface; 
assumption  is  made,  however,  that  flow  within  the  combustion  zone  is  pure 
ial,  i.e.,  zero  axial  component  to  all  ordrs.  Thus  it  could  be  anticipated 
t  the  results  resemble  (regardfng  nonsteady  combustion  behavior)  those  of 
en  [16],  and  there  seems  to  be  only  small  differences  between  the  response  to 
igent  and  to  perpendicular  wave  incidence.  The  problem  is  finally  solved 
lerically,  and  details  of  the  inner/outer  matching  process  were  not  given. 

In  the  remainder  of  this  paper,  the  viscous,  injected  wall  layer 
:mulation  is  derived,  in  perturbation  form.  Analytical  near  field  solutions 
‘  obtained  tot  all  variables  up  to  second  order,  regarding  the  radial 
srdinate  dependence;  the  remaining  (x,t) -dependence  is  shown  to  be  governed  by 
Latively  simple  partial  differential  system.  These  solutions  are  discussed, 
th  particular  attention  to  the  resultant  (first  order)  pressure  distribution 
3  wall  shear  stress,  for  which  experimental  data  are  available. 

AS ALTS IS 


glytical  Model  of  the  Coreflow 

The  equations  of  motion  pertaining  to  the  core-flow  simulation  are 
esented,  for  an  axisymmetric  flow  field.  The  objective  is  to  simulate 
e  cold-flow  experiments  of  Dr.  Brown  at  UTC/CSD,  which  utilir.e  cylindrical 
ometrv.  For  the  coreflow  region,  with  typical  injection  Reynolds  numbers  of 
der  104,  we  assumed  constant  and  uniform  thermophysical  properties.  A 
hematic  of  the  physical  configuration  is  shown  in  Fig.  1.  Turbulence  and 
mbustion  are  precluded  from  the  present  formulation,  for  reasons  discussed 
rlier.  Other  than  these  simplifications,  the  full  compressible,  nonsteady, 
scous  equations  of  motion  are  considered,  with  all  the  dissipative  terms 
eluded . 

The  five  eouations  of  motion,  for  continuity,  radial  momentum,  axial 
imentum  and  energy  are  presented  in  differential  form.  A  caloric  equation  of 
ate  (pertaining  to  perfect  gas)  completes  the  model  to  form  closure  of  the 
pendent  variables. 

The  following  dimensionless  independent  variables  are  introduced,  based  on 
e  two  physical  scales  of  reference  inner  chamber  radius,  R0*,  and  reference 
ijection  velocity,  vQ* ; 

r  ■  r*/R0*,  x  ■  x*/Rq*,  t  -  t*/t0*  (1) 


ere  t0*  -  R0*/v0*  (2) 

ie  dependent  variables  are: 

5-Vfo*.  v  “  v*/v0*,  u  «  U*/V0*,  hQ  «  h*/hc* 
id  p  »  pVPc*  (3 


In  the  last  equations,  the  properties  used  for  non-dimensionalitation 
e  the  reference  (injected)  density,  C*  ,  and  the  reference  chamber  pressure, 
>*;  the  corresponding  thermal  enthalpy,  hD* ,  is  calculated  from  the  caloric 
luation  of  state, 


(4) 


iere  ~t  «Cp/Cv  is  the  specific  heat  ratio,  is  considered  as  constant. 
?ference  speed  of  sound  is 


\J(7-C)hZ 


The 


(5) 
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eing  inviscid,  could  obtain  a  solution  for  the  axial  velocity  component  which 
atisfied  the  no-slip  boundary  condition  at  the  wall.  The  solution  which 
atisfies  the  boundary  data,  namely,  u(x»0)“0,  u(r"l)»0,  and  v(r»l)«l,  yields: 

U  =  itxCosCfr1-) 

if  the  general  family  of  solutions  obtainable,  only  that  which  allows  full 
letermination  of  the  vorticity  (the  azymuthal  component  alone  remains,)  by  the 
ivailable  boundary  data,  is  physically  meaningful;  the  rest  were  therefore 
ejected.  The  axial  pressure  distribution  obtained  from  the  momentum  equation 
s  parabolic, 

[  P„  -  POO  3/f4= 

•his  type  of  injected  flow  field  has  been  investigated  previously  both 
ixperlmentally  and  theoretically.  In  particular  the  early  theoretical  work  of 
lerman  [2],  who  arrived  at  a  power-series  solution  to  the  perturbation  problem 
>f  suction  in  a  flat,  porous-walled  channel,  with  the  suction  Reynolds  number 
lerving  as  small-perturbation  quantity.  The  analytical  results  of  G.  Taylor  [3) 
ind  Wageman  and  Guevara  [4]  more  closely  resemble  the  cosine  terms  of  Culick  11] ; 
joth  [3,4]  have  carried  out  experiments  as  well,  and  both  demonstrated  very  good 
lareement  between  the  measured'  axial  velocity  profiles  and  the  calculated  ones. 

[t  appears  that  Culick  [1]  has  arrived  at  his  results  independently,  since  no 
reference  was  made  to  any  of  the  previous  works.  In  the  experiments  by  Dunlap, 
Willoughby  and  Hermsen  [5],  the  formulation  derived  by  Culick  [1]  was  used  to 
:orrelate  the  measured  data,  again  with  considerable  success,  regarding  the 
roreflow  axial  velocity  profile,  that  is,  away  from  the  close  neighborhood  of 
the  wall. 

Other  experiments  by  Olson  and  Eckert  [6J  and  later  by  Huesman  and  Eckert 
[7]  tend  likewise  to  verify  the  validity  of  this  formulation,  in  particular 
regardina  the  radial  velocity  profile,  which  indeed  exhibits  a  peak  near  the 
sorous  surface  [6],  as  well  as  the  axial  pressure  distribution  (the  latter  shown 
ss  a  linear  correltaion  between  the  friction  coefficient,  Cf,  and  the  inverse 
riean  axial  velocity,  which  are  both  proportional  to  1/x. 

The  recent  (and  ongoing)  experimental  study  by  Brown,  et  al  [8]  provides 
valuable  information  regarding  the  steady  state  axial  pressure  profile  and  the 
sxial  velocity  distribution,  as  well  as  nonsteady  wall  heat  transfer  (obtained 
by  exciting  the  standing  acoustic  modes  in  the  tube).  Departure  of  the  steady 
state  data  from  the  predictions  of  the  aforementioned  formulation  by  Culick  [1] 
was  attributed  to  possible  transition  to  turbulence.  As  will  be  shown  in  this 
study,  the  Dressure  data  obtained  can  be  simulated  very  well  with  a  first-order 
pressure  perturbation,  arising  from  the  laminar  viscous  wall-layer  analysis. 

Earlier,  Yagodkin  [9]  reported  an  experimental  cold  flow  setup,  with  an 
injected  porous  pipe.  The  maximal  injection  Reynolds  number  was  250,  which  is 
2-3  orders  of  magnitude  less  than  that  corresponding  to  actual  internal  rocket 
flows.  Hot-wire  anemometry  was  used  to  obtain  axial  velocity  and  axial  velocity 
fluctuation  vs  axial  and  radial  distance.  Turbulence  intensity  seems  to  peak 
bear  the  surface,  and  decrease  toward  the  centerline  and  toward  the  pipe  wall. 

These  observations  are  Qualitatively  similar  to  those  obtained  later  by  Yamada, 
et  al  [10],  Although  a  transition  region,  at  Reo»100-150,  was  speculated  [9]  to 
involve  "large  eddy  structures",  no  such  evidence  appears  in  the  experimental 
data  reported  [9] . 

Further  studies  by  Yagodkin,  with  Varapaev  [11]  and  Sviridenkov  [12]  are 
theoretical,  and  address  the  problem  of  laminar  stability  of  injected  channel 
flows,  i.e.,  transition  to  turbulence.  Thus,  modified  versions  of  the  Orr- 
Sommerfield  problem  were  investigated  analytically  [11]  and  numerically  [12]. 

Two  related  laminar  flow  stability  analyses  are  by  Goldshtik,  et  al  [13]  and 
Alekseev,  et  al  [14].  None  of  these  theoretical  analyses  indicates  the  presence  > 

of  large  turbulent  eddy  structures  prior  to  a  full  transition  point,  neither  do 
they  obtain  an  origin  of  such  turbulence  on  the  centerline  upstream. 
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FLOID-DYNAMICALLY  COOPLED  SOLID  PROPELLANT 
COMBUSTION  INSTABILITY  -  COLD  PLOW  SIMULATION* 

Dr.  Moshe  Ben-Reuven 

Princeton  Combustion  Research  Laboratories,  Inc. 

Monmouth  Junction,  New  Jersey 

ABSTRACT 

This  analysis  is  aimed  at  the  near-wall  processes  in  an  injected, 
ixisymmetr ic,  viscous  flow.  It  is  a  part  of  an  overall  study  of  solid 
srooellant  rocket  instability,  in  which  cold  flow  simulation  is  evaluated  as  a 
tool  to  elucidate  possible  instability-driving  mechanisms.  One  such  prominent 
nechanism  seems  to  be  visco-acoustic  coupling,  as  indicated  by  earlier  detailed 
srder  of  magnitude  analysis.  The  major  component  of  the  overall  study,  not 
reported  herein,  involves  numerical  simulation  of  the  full  coreflow  equations  of 
notion  (nonsteady,  axisymmetr ic)  by  a  modified  MacCormack  integration  technique. 

The  formulation  is  presented  in  terms  of  a  singular  boundary  layer  problem,  with 
detail  (up  to  second  order)  given  only  to  the  near-wall  region.  The  injection  J 

Revnolds  number  is  assumed  large,  and  its  inverse  square  root  serves  as  an 
appropriate  small-perturbation  quantity.  The  injected  Mach  number  is  also  small, 
and  taken  of  the  same  order  as  the  aforesaid  small  quantity.  The  radial- 
dependence  of  the  inner  solutions  up  to  second  order  is  solved,  in  polynominal 
form.  This  leaves  the  (x,t)  dependence  to  much  simpler  partial  differential 
eouations.  Particular  results  demonstrate  the  existence  of  a  first-order 
pressure  perturbation,  which  arises  due  to  the  dissipative  near-wall  processes. 

This  pressure  and  the  associated  viscous  friction  coefficient  are  shown  to  agree 
very  well  with  experimental  injected  flow  data. 

INTRODOCTON 

This  is  part  of  a  study  aimed  at  elucidation  of  the  physical  mechanisms 
capable  of  driving  acoustic  instability  in  solid  propellant  motors,  particuarly 
of  the  type  termed  velocity-coupled  instability.  Previous  studies  on  the 
coupling  between  velocity  oscillations  and  the  combustion  process  in  solid 
propellant  motors  have  demonstrated  the  complexity  of  the  overall  phenomenon, 
but  have  not  yet  defined  the  basic  mechanisms  nor  how  they  operate  under  flow 
conditions  orevailing  in  rocket  chambers.  Critical  literature  review  and  order 
of  magnitude  analyses  of  velocity  coupling  mechanisms  have  been  carried  out, 
ncludino  visco-acoustic  coupling  and  turbulence  combustion  coupling.  The  major 
goal  of  the  study  is  the  analytical  simulation  of  the  interior  flow  field  within 
a  solid  propellant  qrain.  The  focus  is  on  the  Stokes  layer,  with  the  objective 
of  investigating  the  particular  instability  mechanism  of  visco-acoustic 
coupling.  Preliminary  analysis  has  indicated  that  this  mechanism  is  both 
plausible  and  sufficiently  powerful  to  drive  nonlinear  vibrations;  it  has  been 
shown  that  the  f reouency-dependent  surface  heat  feedback  component,  due  to 
viscous/acoustic  coupling,  has  both  phase  amplitude  ranges  which  would  enable 
driving  of  acoustic  vibrations;  its  amplitude  tends  to  increase  as  the  mean 
coreflow  Mach  number  and  the  frequency  become  higher.  A  comprehensive 
analytical  model  of  the  flow  field  within  the  viscous  wall  layer  region  has  been 
derived,  for  an  axisymmetr ic,  nonsteady  flow  field  configuration.  For  a 
simulation  of  the  cold  flow  test  results  gnerated  at  UTC/CSD,  four  conservation 
equations  are  incorporated  for  continuity,  momentum  and  energy. 

The  major  asoect  of  the  near-wall  behavior  from  the  visco-acoustic  point  of 
view,  is  the  laminar  dissipative  processes  typical  to  that  region.  This 
analysis  is  focused  on  the  near-wall  processes.  Although  the  solultions  derived 
are  nonsteady  in  general,  the  radial  wall-layer  distributions  obtained  could 
best  be  demonstrated  at  steady  state.  For  this  reason,  the  review  herein  is 
limited  to  steady  behavior. 

Culick  (1)  derived  a  solution  to  the  Stokes  stream  function  eouation  or 
flow  in  a  pipe  with  injected  sidewalls.  The  flow  is  rotational,  and  despite 
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FIGURE  3.  DIMENSIONLESS  FIRST  ORDER  AXIAL  PRESSURE 

DISTRIBUTION,  FROM  THE  NEAR-FIELD  ANALYSIS 
HEREIN.  THE  CURRENT  RESULT  IS  PLOTTED  OVER 

the  original  figure  of  csd/utc  (19825,  nhjch 
also  shons  the  parabolic  formula  of  culick. 
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ae  u»2x  was  used,  for  a  cylindrical  port,  and 
subscript  zero  denotes  2eroth  order  convention. 

Now,  front  Bqs.  (29)  and  (62' , 


This  parameter  is  plotted  against  l/2x  (which 
denotes  the  ratio  of  blowing  to  mean  axial 
velocity)  in  Fig.  5.  A  nearly  linear  relationship 
is  obtained,  using  the  coefficient  values  obtained 
from  the  two  data  groups  in  Fig.  3.  In 
comparison,  the  data  obtained  by  Olson  and  Eckert 
[6]  is  considered.  Ref.  6  includes  a  plot  of  the 
ratio  of  (axial  pressure  gradient) /(mean  dynamic 
axial  head)  vs  v0*A^*  -  l/2x.  This  obtains  an 
almost  linear  correlation,  as  would  be  expected 
from  a  parabolic  pressure  drop.  The  slight 
curvature  however,  particularly  apparent  at  small 
values  of  l/2x  <  0.01,  can  be  followed  only  with 
the  present  formulation,  not  with  any  parabolic 
pressure  profile.  Thus,  the  first  order  pressure 
distribution.  Obtained  from  the  viscous  wall  layer 
analvsis,  agrees  well  with  the  measured  data  of 
Brown,  et  al  (8),  while  the  associated  wall 
friction  coefficient  follows  the  same  trend  as 
that  measured  by  Olson  and  Eckert  (6J. 

GCNaUSKMS 

A  derivation  of  the  viscous  wall  layer  regime 
has  been  presented,  pertaining  to  injected  flow  in 
an  axial  porous  tube,  in  simulation  of  interior 
solid  propellant  rocket  flows. 

Solutions  far  the  radial  coordinate  (or  y- 
dependence)  of  all  the  dependent  variables  ip  to 
*  the  first  order  have  been  generated,  in 

polynominal  form.  The  (x,t)-dependence  is  defined 
in  terms  of  a  relatively  simple  partial 
differential  system. 

Particular  results  of  the  analysis  for  the 
special  ca9e  of  steady  state,  are:  (1)  the  first 
order  pressure  perturbation  was  solved  for  uniform 
B^,  and  its  axial  distribution  is  given 
explicitly;  this  term  is  entirely  due  to  the 
laminar  dissipative  wall-layer  processes,  and  (2) 
the  blown  wall  friction  coefficient  was  likewise 
defined.  Both  results  correlate  well  the 
available  e^>erimental  data.  Finally,  (3)  the 
zeroth  order  axial  velocity  distribution  within 
the  layer  Is  linear  radially;  thus,  to  lowest 
order,  viscous  dissipation  is  negligible  in  the 
axial  momentum  balance.  This  indicates  why 
inviscid,  rotational  solutions  (sixh  as  those  of 
Culick  [1]  and  others,  chosen  so  as  to  satisfy  the 
no-slip  condition  at  the  wall)  are  so  successful 
in  representing  this  family  of  flows  -  to  zeroth 
order. 
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The  experimental  data  of  Brown  [8]  also 
demonstrates,  as  evident  in  Fig.  3,  that 
significant  departures  from  the  parabolic  axial 
pressure-drop  profile  (as  predicted  by  the 
inviscid  formulation)  evolve  at  sufficiently  large 
x.  in  the  meantime,  departures  from  the  self¬ 
similar  cosine  velocity  profiles  are  evident, 
which  become  more  appreciable  with  increasing  x. 
Elucidation  of  the  first-order  pressure 
pertrubatian  seems  therefore  highly  important  for 
proper  understanding  of  this  type  of  Injected 
flows.  This  is  undertaken  in  the  remainder  of  the 
present  section,  for  steady  state. 

To  resolve  the  axial  variation  of  p^  by  the 
wall  layer  formulation,  the  second  compatibility 
condition  in  Eg.  (51)  can  he  used,  corresponding 
to  the  y-term: 


DISCUSS  KW  CP  RESULTS 

1t>  facilitate  comparison  with  the 
ejqser imental  data  reported  by  Brown,  et  al  [8], 
one  may  form  the  normalized  axial  pressure  drop, 

)=  MS.'M  *1 

(64) 

This  formula  is  used  to  correlate  the  experimental 
data  of  Brown,  et  al  (8),  as  shown  in  Fig.  3. 
Clearly,  the  measured  pressure  profile  is 
correlated  very  well  by  Eg.  (64),  which  is 
obviously  superior  to  the  inviscid  egression  (1), 
shown  as  well. 


Cj^cCc  ^ 

Fo  TR  -7JK  ■ 


(58) 


For  the  special  case  of  uniform  (zeroth  order) 
injection  at  steady  state,  the  presence  of  a 
nonzero  first-order  pressure  perturbation  would 
imply  physically  a  corresponding  nonzero 
perturbation  upon  the  mass  flux  injected,  i.e.. 


B0  -  F^O.x.t)  /  0 


A  single  point  of  the  data  (x;  p^)  has  been 

utilized  to  obtain  a  scale  for  the  comparison 
(this  is  necessary,  since  no  physical  input  is 
available  regarding  the  value  of  B0,  the 
injected  mass  flux  perturbation,  necessary  for 
defining  bj,  t^),  aJLong  with  pD«l,  F  wvo—l,  and 
K  -  1.4.  suppose  now  that  R^-l,  and  we  select  a 
value  of  B_^0.  (This  is  based  on  some  trial  and 
error  -  but  shows  how  the  correlation  was  obtained 
without  any  regression  analysis);  then. 


as  given  by  Eq.  (42).  Now  at  steady  state, 
although  BQ  is  expected  to  vary  with  d,,  we  will 
(as  a  first  step),  assume  for  simplicity  that 
B0(x)  «  "§0  -  const. 

With  the  foregoing  steady  state  assumptions, 
Bq.  (58)  can  be  readily  turned  into  an  ordinary 
differential  equation,  for  0  <,  x  <  L  s 


where  uniform  injection,  FQ-const,  was  assumed; 
the  coefficients  are: 


bo«l/Bo«l/60,  bj-l/X Bq-1/1 . 4x60 ,  -0.012 
-  1/\/Fbo  -  0.014 

Tt>  demonstrate  the  sensitivity  of  the  axial 
pressure  drcp  to  Reynolds  number  (or  to  the 
parameters  and  &),  the  experimental  points  of 
Fig.  3  were  converted  from  their  original  scale, 

Mi , 

to  the  Ap]/£-fcrm  herein,  using  the  relevant 
values  of  M0  and  R^  given.  The  results  of  this 
conversion,  for  two  distinct  data  groups,  are 
plotted  in  Fig.  4,  along  with  their  excellent 
correlation  by  Bq.  (64) . 


_  b,Ve. 


>0 

j 


nr?a  >  -  (so) 


Note  that  at  steady  state,  according  to  Eqs.  (39) 
and  (44),  Cj-Covo»0.  The  houndary  data  are. 


dpj/dx(0)»O,  and  p^lx-U-p^  (61) 


The  solution  is  straightforward;  for  oositive 
and  bg: 

^(y&b,  *) 


(62) 


fax)  *  b(^|^<^  ^)|(63) 

This  concludes  the  derivation  of  the 
injected,  viscous  wall  layer,  up  to  second  order. 
Full  solutions,  namely,  matcHng  between  inner  and 
outer  expansions  will  not  be  attempted  herein. 
Important  insights  are  obtained  already  from 
resolving  the  near-field  behavior  up  to  second 
order.  In  terms  of  the  y-polynominals. 


In  Fig.  4,  the  differences  between  high 
(low  R^)  and  low  K_  (high  R^,)  are  amplified. 

High  R^_  measurements  appear  on  top,  closest  to 
the  inviscid  pressure  profile. 

Two  important  observations  can  therefore  be 
made:  (1)  axial  pressure  variation  to  lowest 
order  is  0(£),  and  is  governed  by  the  dissipative 
wall  layer  processes,  as  derived  In  the  analysis 
herein.  The  behavior  obtained  in  x  differs  from 
the  parabolic  pressure  drop  formula  of  Culick  (1); 
also,  (2)  one  need  not  invoke  local  turbulence 
generation  to  explain  the  departure  of  measured 
Pl  from  the  predicted  inviscid  behavior. 

Another  property  of  interest  is  the  wall 
friction  coefficient,  or  dimensionless  wall  shear 
stress, 

Cf  =  &/ 

where  u*  denotes  the  mean  axial  ooreflow  velocity. 
Using  dimensionless  convention  employed  herein, 
along  with  the  wall  layer  coordinate, 

Cfo  =  fc  (2Y*L  (65) 
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Thus,  h^(y,x,t)  is  second  order  in  y,  like 
the  appropriate  boundary  condition  is, 

B2<x,t)  -  hj^O.x.t)  -  h^Bo  (48) 

The  axial  momentum  balance,  Bq.  (33),  after 
dividing  through  by  FQ  and  following  integration, 
yields: 


and 


(49) 
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undetermined  coefficients  should  arise 
necessarily,  to  accommodate  coupling  with  the 
outer,  irrviscid  (core)  flowfield,  through  inner/ 
outer  asymptotic  matching. 


THE  FIRST  OREER  PPES5UBE  PERTURBATION 

In  the  outer  inviscid  flow  (coreflow),  the 
pressure  difference  is  balanced  tv  the  axial 
acceleration,  as 


Uj^O.x.t)  -  0 

satisfying  the  no-slip  condition  at  the  wall. 

Thus,  the  perturbed  axial  velocity  is  third  order 
in  its  y-dependenoe,  and  the  corresponding  viscous 
dissipation  term  (unlike  its  zerotb-order 
counterpart),  does  not  vanish  within  the  layer. 

With  the  foregoing  derived  results,  the 
zeroth  order  energy  equation  can  be  shown  to  yield 
merely  a  condition  connecting  UQ  and  vQ  timewise. 
One  mav  turn  now  to  the  first-order  energy 
equation,  which  seems  to  yield  same  simple  and 
highly  useful  results  even  without  full  solution. 
After  some  manipulation,  Bq.  (38)  obtains: 


In  the  meantime,  the  available  CSD  experimental 
data  clearly  shows  that  the  seemingly  equivalent 
representation, 

<X.h)  ,M, 

where  0  <  £«  1.  Therefore,  within  the  range  of 
flow  conditions  simulated  herein,  it  is  inferred 
that  the  mean  axial  coreflcw  velocity  can  become 
considerably  larger  than  the  radial  (or  injected) 
velocity,  such  that,  from  Eqs.  (52)  and  (53): 

U.*(OUTfK )/)/*(  (Hj£CT»0^)  /v  0(_  m 

and  furthermore. 


-  "X=i .  ss  o 

■znc  rT  2x3*- 


(50) 


*  The  first  bracketed  term,  after  using  the 
foregoing  results,  is  simply 


(P* (0)-p*(x) ]/p* (x)  "  ~  0(£). 

Recall  that  pQ  «  pQ(t)  and  p-,  «  p^Xft)  as 
derived  in  wadi  layer  analysis.  Hence,  up  to 
first  order,  no  radial  pressure  gradients  arise, 
and  the  inner  and  outer  pressures  must  be  equal, 
viz., 


(o) 


(i) 


Pi 


(o) 


Pi 


(i) 


(55) 


vo  -  ^r^o'^o 

Using  the  appropriate  first  order  expressions 
obtained  herein  (for  u^,  v^,  and  hjj  in  Eq.  (50), 
and  collection  of  equal  powers  of  y  yields: 


+ S<v«'  SP )- 

Yp,{'  Jt(  Ml  l  y  + 

°  F<?  3 


(51) 

Oompatibility  with  the  foregoing  derivation 
(In  which  y  and  (x,t)  variable  separation  was 
Implemented),  can  be  maintained,  provided  each  of 
the  bracketed  berms  in  Bq.  (51)  vanishes 
identically.  The  resulting  four  oompatibility 
relations, in  psrtial  differential  form, would 
determine  the  behavior  of  the  wall  suhlayer  system 
up  to  the  first  order  In  &  ,  the  small 
perturbation  quantity.  However,  a  total  of  three 


Suppose  now  that  the  leading  term  in  the  outer 
axial  velocity  expansion  is,  according  to  Bq. 

(54), 

(KI 

while  the  remaining  outer  variables  are  of  simpler 
form  $  e «g • f 

v(o)(x/,t)  ~  V0(0)+  V,(<V  Kttt. 

Therefore,  the  outer  axial  momentum  balance, 
derived  from  Eqs.  (7)  —  (12)  yields,  at  order  1/fc. : 

which  shows  that  indeed  the  outer  flow  can  support 
this  first-order  pressure  perturbation.  To 
conclude,  a  first-order  pressure  perturbation 
within  the  flow  field,  £,p,(x,t),  has  been 
postulated,  following  the  viscous  wall  layer 
analysis.  This  pressure  is  common  to  both  inner 
(wall  layer)  and  outer  (coreflow)  regions.  Within 
the  wall  layer,  it  is  balanced  by  the  shear  force, 
or  zeroth-order  vocticity  generation,  as  shown  in 
Eq.  (28).  In  the  coreflow,  it  is  balanced  by  the 
lowest-ocder  axial  acceleration. 
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4* 


poeltion  (x,t).  Even  more  striking  is  the 
vanishing  of  the  viscous  dissipation  term  at 
zeroth  order: 


ss  o 


whidi  leaves  in  the  zeroth-order  axial  momentum 
equation  a  balance  of  inertial  terms,  strictly. 
This  explains  physically  the  success  (up  to  first 
order)  of  modeling  this  family  of  injected  flows 
by  assuming  rotational,  inviscid  motions:  such 
modeling  indeed  obtains  solutions  for  the  axial 
velocity  profile,  which  satisfy  the  no-slip 
condition  at  the  wall  (r«l). 

It  further  appears  that  the  shear  stress,  Bq. 

(30),  is  proportional  to  the  first-order  axial 
pressure  gradient,  while  being  inversely 
proportional  to  the  injected  mass  flux,  as  would 
be  expected.  Of  course,  depends  on  F_,  and 

one  expects  their  ratio  to  be  finite  at  the  limit 
as  zero  injection  is  approached. 

The  zeroth-order  differential  system  reads: 


aft 


~D  _ 

*5x 


- -  tr 

—  *  a 


Or  *  *  -- 


(31) 


(32) 


(33) 


+^C^^V^C-fF;v0^vO=  -^v0+ 

tr-o 


reads: 


The  corresponding  first-order  formulation 


(34) 


(35) 


'3’fi  _ rr  _  t,  r- 

(36) 

+fiv4ui)*" 

(37) 


H  ^  V,  )=  -tPM  + 

|k  )  +- 


(38) 


With  the  foregoing  result  for  axial  velocity, 
the  zeroth  order  formulation  can  be  utilized  to 
solve  for  the  y-dependence  of  the  other  dependent 
variables.  The  zeroth  order  continuity  equation 
can  be  written  in  split  form,  since  Fc  and  are 
independent  of  y: 


'zz/dfc  +  F.  -  c0(*-jt) 

2.  (*&  VB.  V.  _  IfF,  _  c  .  , 

v.  tix™  ^vj  ~  ~  C.(*>v 


(39) 


(40) 


where  CQ(x,t)  is  a  common  separation  parmeter, 
with  a  range  of  values  uniquely  corresponding  to 
the  boundary  data.  The  second  equation  yields 

Fj^x.t)  «  ^(x.t)  +  CQ(x,t)  y+^flfejUOy2  (41) 

where: 

B0(x,t)  «  Fj^O.x,^ 
u0(x,t>  s 

'  ho 

Similarly,  the  zeroth  order  radial  momentum 
equation  yields,  after  splitting: 

Qf*(dt+  =  C,o,*) 


(42) 

(43) 

(44) 


(45) 


JL-ir .  -L-av.  J* 

ei(x;t)=  v^tjO) 

The  foregoing  results  for  Fj  and  Vj^  yield  for  the 
first-order  density:  A 

^y)=^o-fe6,yv0+  (2c«-c,/v0v 

(46) 

Note  that:  fx(0,x,t)  -  ©0  -  f^v,, 

so  that  is  uniquely  defined  and  an  additional 
integration  constant  is  not  necessary.  Now, 
following  the  definition  of  and  since 

=0, 

then: 

hj (y»x,t)  -  -^^>(Y>x,t)  +  BjfXjt)  (47) 
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<c)  At  the  (norpermeable  solid)  head-end  closure, 
(t,  r,  x»0): 


v  -  0, 


0,  h»hH(r,t) 


The  functions  vQ(x,t),  )u(x,t)  and  fy(r,t)  are 
arbitrary  imposed  distributions. 

(d)  The  exit  plane,  defined  by  (t,  r,  x«L),  forms 
an  entrance  into  a  short,  convergent  nozzle 
section. 

THE  mJBCTED  SIDEWALL  LAYER 

The  flow  region  of  interest  is  close  to  the 
surface,  where  viscous  forces  are  expected  to  be 
appreciable  within  a  thin  layer. 

For  the  neighborhood  of  r«l,  the  following 
transform  is  proposed  for  the  radial  coordinate: 

y  -  (l-r)/£  (17) 

which  magnifies  the  wall  layer,  with 


0  <  £  *  1//^  «  1 


y>r-kb> 


where  is  the  relevant  (injection)  Strouhal 
number.  The  range  of  Strouhal  numbers  considered 
is  "  0(1),  90  one  need  not  introduce  any 
additional  timescales. 

The  independent  variables  are  (x,y,t),  while  the 
associated  dependent  variables,  in  the  wall-layer, 
are  perturbed, 

,  v=v0-i-fev0  u  =  T^o+-e.w»; 

k0+fek,  ,  (24) 

while  the  following  abbreviations  are  introduced, 

Fe  «  To V.  ,  F,  S  +?,  \Jo 

Sc=?,u«  ,  G,s  <25) 

It  should  be  stresed  that  the  0(  £, )  berms 
represent  perturbation  quantities  which  may  be 
later  considered  as  series  expansions.  The 
present  analysis  is  concerned  with  the  two  lowest 
orders  only.  This  in  no  way  implies  limitation 
to  so-called  "linear"  considerations. 

Far  the  perturbation  variables  of  Egs.  (17)- 
(25),  the  continuity  equation  yields 

4(g.tES,)  =  -(l-6y)(Vf<0 


The  assumption  for  small  injection  Mach  number  is 
constrained  as  follows.  Obviously,  the  injection 
Mach  number  appears  as  an  additional  parameter  in 
the  formulation  (equations  of  momentum  and 
energy).  In  the  flow  types  of  interest  for 
simulation  herein,  Mp  is  also  very  small;  in 
consideration  of  typical  experimental 
configurations  at  CSD/UTC,  we  find 


^-oO/RecVOCO 


which  adequately  represents  a  range  of  cold- flew 
conditions.  This  offers  great  simplification  in 
the  analysis,  although  at  the  cost  of  narrower 
range  of  general  application  (considering  the 
relative  freedom  of  the  two  major  flow  parameters, 
Rgp  and  M^).  Therefore,  a  parameter  is 
introduced. 


__  vjveo  K - a/,1 

-  DM*  _  U{  ’  <20 

The  question  of  timescale  depends  upon  the 
range  of  frequencies  of  interest.  The  following 
reasoning  will  demonstrate  that  for  the  range  of 
conditions  considered  for  the  present  simulation, 
the  timescale  can  remain  the  same  one  as  in  the 
ooreflow.  The  viscous  layer  thickness  is 

*v  -  v/^T  <21 

The  Stokes  Layer  thickness  (for  acoustic 
perturbations  with  a  frequency  f0) 

*310  “(?/ <22 

The  ratio  of  these  two  thickness  scales  is, 

Sy  6sro  ‘/  "w/v*-  Jho  <23 


Similes  substitution  of  perturbed  variables  is 
also  carried  out  for  the  remaining  equations  of 
motion;  a  detailed  derivation  is  available  in  Ref. 
17.  A  hierarchy  of  equations  can  then  be 
collected,  for  equal  powers  of  the  small  quantity^. 

The  lower-order  analysis  (concerning  negative 
powers  of  £ )  readily  yields  the  followinq  simple 
results: 

P0  "  F0(x,t)  *  vo  “  vo<x'fc>  '  Po  “  Po*^ ' 

Pi  “  Pi(x,t)  (27) 

Also,  the  following  differential  equation  arises 
from  the  axial  mcmentisn  balance  at  order  1/  £  : 


which  can  be  integrated  to  reveal  the  y-deperrience 
of  Ug: 

with  (x,t>0)  ■  0  . 

Of  course,  the  (x,t)  dependence  of  Uq  still 
remains  to  be  four*!  However,  Its  dependence 
upon  the  layer  coordinate,  y,  is  found  to  be 
linear;  tills  result  has  several  important 
implications.  The  axial  shear  stress  component 
within  the  layer, 

_.o  ^  _  K* 

~  ^  f=„  7PC  ,  (30) 

is  obviously  nonzero  in  general,  while  being 
Independent  of  distance  from  the  wall  at  any  given 
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Analytical  Model  of  the  Coreflow 

IT*  aquations  of  motion  pertaining  to  the 
core-flow  simulation  are  presented,  for  an 
axlsymmetric  flow  field.  The  objective  is  bo 
simulate  the  cold- flow  experiments  of  Dr.  Brown  at 
tmyCSD,  which  utilize  cylindrical  geometry.  For 
the  coreflow  region,  with  typical  injection 
Reynolds  numbers  of  order  10®,  we  assumed  constant 
and  uniform  thermophysical  properties.  A 
schematic  of  the  physical  configuration  is  shown 
in  Fig.  1.  Turbulence  and  combustion  are 
precluded  from  the  present  formulation,  for 
reasons  discussed  earlier.  Other  than  these 
simplifications,  the  full  compressible,  nonsteady, 
viscous  equations  of  motion  are  considered,  with 
all  the  dissipative  terms  Included. 


These  idealizations  are  incorporated  merely  for 
convenience,  in  allowing  clear  identification  of 
physical  interactions  within  the  ooreflow,  at  low 
(axial)  Mach  numbers;  sharp  pressure  and 
temperature  variations  are  otwiously  precluded. 
The  dimensionless  equations  of  motion  are  as 
follows,  for  the  region 

0  <  x  <  1,  0<r<l,  t>0: 


■frtr 


(7) 


where  the  dependent  variable  vector  is, 

=  , *v>  (8) 
The  radial  and  axial  flux  terms  are,  respectively. 


The  five  equations  of  motion,  for  continuity, 
radial  momentum,  axial  momentum  and  energy  are 
presented  in  differential  form.  A  caloric 
equation  of  state  (pertaining  to  perfect  gas) 
completes  the  model  bo  form  closure  of  the 
dependent  variables. 


pT  ■  ?«*v,  {9) 

gt-  (fu,  ,rph<0 

(10) 


The  following  dimensionless  independent 
variables  are  Introduced,  based  on  the  two 
physical  scales  of  reference  inner  chamber  radius, 
Rjj*,  and  reference  injection  velocity,  vQ*; 

r  -  r*/Rc*,  x  «  x*/R0*.  t  -  t*/tQ*  (1) 
where  V  -  R^/V 


where  superscript  T  denotes  vector  transpose.  The 
right-hand  side  (source)  terms  are  defined:  Sj-0, 


V}  X- 

tZeo  r  v-wr  r  /  -gr  ~ 


irr’S  Eeo 


(11) 


The  dependent  variables  are: 

v  *  v*/vQ*,  u  «  u*/v0*. 

?=?*/<?*,  k=kVkt, 


The  properties  used  for  non-dimensionalization  are 
the  reference  (injected)  density,  ,  and  the 
reference  chamber  pressure,  PQ*;  the  corresponding 
thermal  enthalpy,  hQ*,  is  calculated  from  the 
caloric  equation  of  state. 


f*  =  *£±  P*ht 

■o  Y  »o 


(3) 


where  «Cp/Cv  is  the  specific  heat  ratio,  is 
considered  a  constant.  The  reference  speed  of 
sound  is 

The  corresponding  injection  Mach  nurfcer  is 


(5) 


The  reference  (injection)  Reynolds  number,  and 
Prandtl  number  are,  respectively. 


£<0  ~ 


r.v«? 


T?  =|U*<VW 


(6) 


Recall  that  the  viscosity,  thermal  conductivity 
and  isobaric  specific  heat  are  all  uniform  and 
constant  within  the  present  oold-flcw  simulation. 


(13) 


The  following  physical  boundary  data  are 
available,  for  the  cold- flow  simulation: 

(a)  On  the  centerline,  (t,  r^),  x) : 

'**>>  11,1 


(b)  At  the  porous  (injected)  surface,  (t,  r-1,  x) 
v  ■  -vc(x,t),  u  •  0,  h  ■  !\,(x,t)  (15) 


The  major  aspect  of  the  near-wall  behavior 
from  the  visoo- acoustic  point  of  view,  is  the 
laminar  dissipative  processes  typical  to  that 
region.  This  analysis  is  focused  on  the  near-wall 
processes.  Although  the  solultions  derived  are 
nonsteady  in  general,  the  radial  wall-layer 
distributions  obtained  could  best  be  demonstrated 
at  steady  state.  For  this  reason,  the  review 
herein  is  limited  to  steady  behavior. 

Culick  [1]  derived  a  solution  bo  the  Stokes 
stream  function  equation  or  flow  in  a  pipe  with 
injected  sidewalls.  The  flow  is  rotational,  and 
despite  beinq  inviscid,  could  obtain  a  solution 
for  the  axial  velocity  component  which  satisfied 
the  no-slip  boundary  condition  at  the  wall.  The 
solution  which  satisfies  the  boundary  data, 
namely,  u(x«0)“0,  u(r«l)”0,  and  v(r«l)«l,  yields: 

V=-SiK£rl)/r 

tL=  nxoo(§-rl) 

Of  the  general  family  of  solutions  obtainable, 
only  that  which  allows  full  determination  of  the 
vorticity  (the  azymuthal  component  alone  remains,) 
by  the  available  boundary  data,  is  physically 
meaningful;  the  rest  were  therefore  rejected.  The 
axial  pressure  distribution  obtained  from  the 
momentum  equation  is  parabolic, 

[f 

This  type  of  injected  flow  field  has  been 
investigated  oreviously  both  experimentally  and 
theoretically.  In  particular  the  early 
theoretical  work  of  Berman  [2],  who  arrived  at  a 
power-series  solution  to  the  perturbation  problem 
of  suction  in  a  flat,  porous-walled  channel,  with 
the  suction  Reynolds  number  serving  as  small- 
perturbation  quantity.  The  analytical  results  of 
G.  Taylor  [3]  and  Wageman  and  Guevara  [4]  more 
closely  resemble  the  cosine  terms  of  Culick  (1); 
both  [3,4]  have  carried  out  experiments  as  well, 
and  both  demonstrated  very  good  agreement  between 
the  measured  axial  velocity  profiles  and  the 
calculated  ones.  It  appears  that  Culick  [1]  has 
arrived  at  his  results  independently,  since  no 
reference  was  made  to  any  of  the  previous  works. 

In  the  experiments  by  Dunlap,  Willoughby  and 
Hermsen  [5],  the  formulation  derived  ty  Culick  [1] 
was  U3ed  to  correlate  the  measured  data,  again 
with  considerable  success,  regarding  the  core flow 
axial  velocity  profile,  that  is,  away  from  the 
close  neighborhood  of  the  wall. 

Other  experiments  by  Olson  and  Eckert  [6]  and 
later  ty  Huesman  and  Eckert  [7]  tend  likewise  to 
verify  the  validity  of  this  formulation,  in 
particular  regarding  the  radial  velocity  profile, 
which  indeed  exhibits  a  peak  near  the  porous 
surface  [6],  as  well  as  the  axial  pressure 
distribution  (the  latter  shown  as  a  linear 
correltaion  between  the  friction  coefficient,  Cf, 
and  the  inverse  mean  axial  velocity,  which  are 
both  proportional  to  1/x, 

The  recent  (and  ongoing)  experimental  study 
by  Brown,  et  si  [8]  provides  valuable  Information 
regarding  the  steady  state  axial  pressure  profile 
and  the  axial  velocity  distribution,  as  well  as 
nonsteady  wall  heat  transfer  (obtained  by  exciting 
the  standing  acoustic  modes  in  the  tube). 


Departure  of  the  steady  state  data  from  the 
predictions  of  the  aforementioned  formulation  by 
Culick  (1)  was  attributed  to  possible  transition 
to  turbulence.  As  will  be  shown  in  this  study, 
the  pressure  data  obtained  can  be  simulated  very 
well  with  a  first-order  pressure  perturbation, 
arising  from  the  laminar  viscous  wall-layer 
analysis. 

Earlier,  Yagodkin  [9)  reported  an 
experimental  cold  flow  setup,  with  an  injected 
porous  pipe.  The  maximal  injection  Reynolds 
number  was  250,  which  is  2-3  orders  of  magnitude 
less  than  that  corresponding  to  actual  internal 
rocket  flows.  Hot-wire  anemometry  was  used  to 
obtain  axial  velocity  and  axial  velocity 
fluctuation  vs  axial  and  radial  distance. 

Turbulence  intensity  seems  to  peak  near  the 
surface,  and  decrease  toward  the  center line  and 
toward  the  pipe  wall.  These  observations  are 
qualitatively  similar  to  those  obtained  later  by 
Yamada,  et  al  [10].  Although  a  transition  region, 
at  RgQ-100-150,  was  speculated  [9]  to  involve 
"large  eddy  structures",  no  such  evidence  appears 
in  the  experimental  data  reported  [9]. 

Further  studies  by  Yagodkin,  with  Varapaev 
[11]  and  Sviridenkov  [12]  are  theoretical,  and 
address  the  problem  of  laminar  stability  of 
injected  channel  flows,  i.e.,  transition  to 
turbulence.  Thus,  modified  versions  of  the  Orr- 
Sommerfield  problem  were  investigated  analytically 
[11]  and  numerically  [12].  Two  related  laminar 
flow  stability  analyses  are  by  Goldshtik,  et  al 
[13]  and  Alekseev,  et  al  [14],  None  of  these 
theoretical  analyses  indicates  the  presence  of 
large  turbulent  eddy  structures  prior  to  a  full 
transition  point,  neither  do  they  obtain  an  origin 
of  such  turbulence  on  the  centerline  upstream. 

Recently  Flan dro  [15]  has  carried  out  a 
theoretical  analysis  for  a  burning  propellant  in  a 
cylindrical  grain,  under  the  effect  of  incident 
acoustic  waves.  A  detailed  formulation  was 
derived  with  a  double  expansion,  in  terms  of  both 
inverse  Reynolds  number  as  well  as  Mach  number 
(independent  small  parameters).  A  nonsteady 
premixed  combustion  zone  was  considered  near  the 
propellant  surface;  the  assumption  is  made, 
however,  that  flow  within  the  combustion  zone  is 
pure  radial,  i.e.,  zero  axial  comonent  to  all 
ordrs.  Thus  it  could  be  anticipated  that  the 
results  resemble  (regarding  nonsteady  combustion 
behavior)  those  of  T'ien  [16],  and  there  seems  to 
be  only  small  differences  between  the  response  to 
tangent  and  to  perpendicular  wave  incedence.  The 
problem  is  finally  solved  numerically,  and  details 
of  the  inner/outer  matching  process  were  not  given. 

In  the  remainder  of  this  paper,  the  viscous, 
injected  wall  layer  formulation  is  derived,  in 
perturbation  form.  Analytical  near  field 
solutions  are  obtained  foe  all  variables  up  to 
second  order,  regarding  the  radial  coordinate 
dependence;  the  remaining  (x,t) -dependence  is 
shown  to  be  governed  by  relatively  simple  partial 
differential  system.  These  solutions  are 
discussed,  with  particular  attention  to  the 
resultant  (first  order)  pressure  distribution  and 
wall  shear  stress,  foe  which  experimental  data  are 
available. 
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ABBIMCr 

This  analysis  is  aimed  at  the  near-wall 
processes  in  an  injected,  axisymmetric,  viscous 
flow.  It  is  a  part  of  an  werall  study  of  90lid 
propellant  rocket  instability,  in  which  cold  flow 
simulation  is  evaluated  as  a  tool  to  elucidate 
possible  instability-driving  mechanisms.  One  such 
prominent  mechanism  seems  to  be  visco-acoustic 
oouoling,  as  indicated  by  earlier  detailed  order 
of  magnitude  analysis.  The  major  component  of  the 
overall  study,  not  reported  herein,  involves 
numerical  simulation  of  the  full  coreflow 
equations  of  motion  (nonsteady,  axisummetrlc)  by  a 
modified  MacCormack  integration  technique.  The 
formulation  is  presented  in  terms  of  a  singular 
boundary  layer  problem,  with  detail  (up  to  second 
order)  given  only  to  the  near-wall  region.  The 
injection  Reynolds  number  is  assumed  large,  and 
its  inverse  square  root  serves  as  an  appropriate 
small-perturbation  quantity.  The  injected  Mach 
number  is  also  small,  and  taken  of  the  same  order 
as  the  aforesaid  small  quantity.  The  radial- 
dependence  of  the  inner  solutions  up  to  second 
order  is  solved,  in  polynominal  form.  This  leaves 
the  (x,t)  dependence  to  much  simpler  partial 
differential  equations.  Particular  results 
demonstrate  the  existence  of  a  first-order 
pressure  perturbation,  which  arises  due  to  the 
dissipative  near-wall  processes.  This  pressure 
and  the  associated  viscous  friction  coefficient 
are  shown  to  agree  very  well  with  experimental 
injected  flow  data. 

NDCNdAIURE 

•  nr'zle  throat  area  and  port  exit  area, 
respectively 

«  adiabatic  velocity  of  sound 

■  wall  friction  coefficient, 

«  isochoric  and  isobar ic  specific  heats 

(JAg-K) 

■  radial  mass  flux  (dimensionless) 

■  axial  mass  flux  (dimensionless) 

•  thermal  enthalpy,  dimensionless 
«  ratio  of  inverse  Reynolds  number  and 

Mach  number  squared 

■  chamber  length 

■  Mach  number 

■  pressure 

■  Prandtl  number 

■  channel  radius 

■  injected  Reynolds  nunber 

■  radial  coordinate 

■  “aource'-terms  in  the  eolations  of 
motion  for  coreflow 

•  Strouhal  number,  injected 
»  time  (dimensionless) 

■  parameter  defining  (x,t)  -  variation 
of  wall  layer  axial  velocity  component 

■  axial  velocity,  and  mean  axial 
ooreflow  velocity  respectively 

^beaearch  Scientist 
Member,  AIM 


v  »  radial  velocity  oomponent 

x  »  axial  distance  (dimensionless) 

y  «  radial,  magnified  wall  layer 

coordinate,  perpendicular  to  surface, 


Greek  Symbols: 

Tf  ■  Cr/Cy  specific  heat  ratio 

«  difference,  increment 

*  length  scales 

&  ■  small  perturbation  quantity 

*  thermal  conductivity  of  gas  (air), 
JA-m-s 

fX  »  viscosity  ooeff icient,  kg/m-s 

<^>  ■  density 

Subscripts,  Superscripts: 

(  )Q  «  denotes  zeroth  order  (perturbation) 

{  )i  «  denotes  first  order  perturbation 

(  )*  ■  denotes  dimensional  quantity 

BACRGBOOND 

This  is  part  of  a  study  aimed  at  elucidation 
Of  the  physical  mechanisms  capable  of  driving 
aooustic  instability  in  solid  propellant  motors, 
particularly  of  the  type  termed  velocitv-ooupled 
instability.  Previous  studies  an  the  coupling 
between  velocity  oscillations  and  the  combustion 
process  in  solid  propellant  motors  have 
demonstrated  the  (complexity  of  the  overall 
phenomenon,  but  have  not  yet  defined  the  basic 
mechanisms  nor  how  they  operate  under  flow 
conditions  prevailing  in  rocket  chambers. 

Critical  literature  review  and  order  of  magnitude 
analyses  of  velocity  coupling  mechanisms  have  been 
carried  out,  including  visoo-acoustic  coupling  and 
turbulence  combustion  coupling.  The  major  goal  of 
the  study  is  the  analytical  simulation  of  the 
interior  flex  field  within  a  solid  propellant 
grain.  The  focus  is  on  the  Stokes  layer,  with  the 
objective  of  investigating  the  particular 
instability  mechanism  of  visoo-acoustic  coupling. 
Preliminary  analysis  has  indicated  that  this 
mechanism  is  both  plausible  and  sufficiently 
powerful  to  drive  nonlinear  vibrations;  it  has 
been  shown  that  the  frequency-dependent  surface 
heat  feedback  oomponent,  due  to  viscous/aooustic 
coupling,  has  both  phase  and  amplitude  ranges 
which  would  enable  driving  of  acoustic  vibrations; 
its  amplitude  tends  to  increase  as  the  mean 
coreflow  Mach  number  and  the  frequency  become 
higher.  A  comprehensive  analytical  model  of  the 
flow  field  within  the  viscous  wall  layer  region 
has  been  derived,  for  an  axisymmetric,  nonsteady 
flow  field  configuration.  For  a  simulation  of  the 
cold  flow  test  results  gnerated  at  UTC/CSD,  four 
conservation  equations  are  incorporated  for 
continuity,  momentum  and  energy. 
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The  corresponding  injection  Mach  number  is 

=  v.Va?  ,6) 

The  reference  (injection)  Reynolds  number  and  Prandtl  number  are,  respectively, 

9.-v.-e:/ju-  ,  -pr=fj.-4 W 

Recall  that  the  viscosity,  thermal  conductivity  and  isobaric  specific  heat  are 
all  uniform  and  constant  within  the  present  cold-flow  simulation.  These 
idealizations  are  incorporated  merely  for  convenience,  in  allowing  clear 
identif ication  of  physical  interactions  within  the  coreflow,  at  low  (axial)  Mach 
numbers;  sharp  pressure  and  temperature  variation  are  obviously  precluded. 

The  dimensionless  equations  of  motion  are  as  follows,  for  the  region 

0  <  x  <  1,  0<r<l,  t  >  0 :  __ 


AA  i- (rF  }  +  _  < 

where  the  dependent  variable  vector  is, 

“T-C?  fkj 

The  radial  and  axial  flux  terms  are,  respectively, 

rT«  f v z  ,  fw  Tfhv  } 


-T  _  / 


(8) 


(?) 


(10) 


(11) 


where  superscript  T  denotes  vector  transpose.  The  right-hand  side  (source)  terms 
are  defined:  S-jsO, 


‘/r 


(12) 


(13) 


(14) 


The  parameters,  Y”,  Pr,  M  ^ ,  R  are  all  constants.  Note  that  incorporation  of 
t“.»  Dressure  gradient  within  tne  R,  source  term  in  the  radial  momentum  equation, 
“■•.le  the  axial  pressure  gradient  is  included  within  the  axial  flux  comonent, 

' i ,  is  merely  for  convenience.  In  the  meantime,  the  viscous  and  thermal 
4 ) • « ; r s * l ve  terms,  ~0(l/Reo),  are  expected  to  be  very  small  over  most  of  the 
■ » ' i -w  domain,  excluding  the  neighborhood  of  the  walls. 


UNCLASSIFIED 


-77- 


UNCLASSIFIED 


The  following  physical  boundary  data  are  available,  for  the  cold-flow 
simulation: 

(a)  On  the  centerline,  (t,  r«0,  x) : 


V  =  0 ?  Bh,far=c> 


(b)  At  the  porous  (injected)  surface,  (t,  r»l,  x) : 

v«-v0(x,t),  u-0,  h  -  hD(x,t)  (16) 

(c)  At  the  (nonpermeable  solid)  head-end  closure,  (t,  r,  x-0): 

v  -  0,  u-0,  h  -  hH(r,t)  (17) 

The  functions  ve(x,t),  hQ(x,t)  and  hH(r,t)  are  arbitrary  imposed  distributions. 

(d)  The  exit  plane,  defined  by  (t,  r,  x-L),  forms  an  entrance  into  a  short, 
convergent  nozzle  section.  This  nozzle  section  is  treated  separately  from 
the  rest  of  the  flow  field. 

The  foregoing  discussion  has  summarized  the  coreflow  analytical  model, 
including  the  equations  of  motion  and  the  relevant  boundary  data." 

THE  INJECTED  SIDEWALL  LAYER 

The  flow  region  of  interset  is  close  to  the  surface,  where  viscous  forces 
are  expected  to  be  appreciable  within  a  thin  layer,  as  shown  i En.  2 

For  the  neighborhood  of  r-1,  the  following  transform  is  proposed  for  the 
radial  coordinate: 


y  -  (l-r)/£  (18) 


which  magnifies  the  wall  layer,  with 

0  <  fc  -  l/f%0  <<  1 


>  ■v'/sr'*  t* 

and  r  -  1  -  £y 

The  assumption  for  small  injection  Mach  number  is  constrained  as  follows. 
Obviously,  the  injection  Mach  number  appears  as  an  additional  parameter  in  the 
formulation  (equations  of  momentum  and  energy).  In  the  flow  fields  of  interest 
for  simulation  herein,  M0  is  also  very  small;  in  consideration  of  typical 
experiments  at  CSD/UTC  with  air  injection,  we  find 

mq2  -  O^/^eo') 

which  adequately  represents  a  range  of  cold-flow  conditions.  This  offers  great 
simplification  in  the  analysis,  although  at  the  cost  of  narrower  range  of 
general  apolication  (considering  the  relative  freedom  of  the  two  major  flow 
parameters,  Reo  and  Mc) .  Therefore,  a  parameter  is  introduced, 

■£«.=  (* =-^C  ~ 0(0 


The  question  of  timescale  depends  upon  the  range  of  frequencies  of 
interest.  The  following  reasoning  will  demonstrate  that  for  the  range  of 


(19) 

?yv  <20, 
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resent  simulat 
he  viscous  lay 

Sv  -  R 0*/V^ 


conditions  considered  for  the  present  simulation,  the  timescale  can  remain  the 
same  one  as  in  the  coreflow.  The  viscous  layer  thickness  is 


'  eo 


(22) 


The  Stokes  Layer  thickness  (for  acoustic  perturbations  with  a  frequency  f0) 

$STO  -C^o)1'3-  (23) 


The  ratio  of  these  two  thickness  scales  is, 

W£sto  *  \ff0V/^“  'Ts 


RO 


(24) 


where  SR0  is  the  relevant  (injection)  Strouhal  number.  The  range  of  Strouhal 
numbers  considered  is  SRQ  “  0(1),  so  one  need  not  introduce  any  additional 
timescales. 

The  independent  variables  are  (x,y,t),  while  the  associated  dependent 
variables,  in  the  wall-layer,  are  perturbed, 

f-fo+Zf,  J  V*V0-t-£Vt  5  u  =  + 


(25) 


while  the  following  abbreviations  are  introduced, 

D  &o-foU-C  ,  ,  ^1=  ,26) 

For  the  perturbation  variables  of  Eqs.  (18)-(26),  the  continuity  equation 
yields 


(27) 


The  following  hierarchy  is  collected: 
ORDER  l/£ 

ORDFP  P°  (ZEROTH) 


,  'c)  <=>© 

KE. 

2Tb 

(FIRST) 

HSi 

"$<2*  - 

'ZPc  ' 

r  zx. 

(28) 


(29) 


(30) 


The  radial  momentum  balance,  after  similar  substitution,  yields  the  following 
hierarchy. 

ORDER  1/  £3: 


0  - — ^  Vo  ~ 


yRDEP  1/  6J-. 


(31) 


(32) 
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ORDER  1/g.  : 


-7>(F0V«)/^  =0 


(33) 


ORDER  fe,°  (ZEROTH) 


H-1-  £  (Few.)-  ^(F.  V,  *F,W)-  -RWt  ,341 


ORDE R _ £ 1  (F I RST)  : 


+  ^v'  +p-<0  . 

Similarly,  the  axial  momentum  balance  yields  the  associated  hierarchy: 

ORDER  1/  E 2 : 


(35) 


ORDER  1/  E  : 


!• 


KrSZPolHoc^  0  =>  pe~p.  (tr) 


K.1&-  =  ° 


(36) 


(37) 


where  we  used  F0u0  •  G0v0.  This  equation  is  of  great  importance,  as  will  be 
shown  later. 


ORDER  E°  ( ZEROTH  1: 

+^0^^  &<^-  'f  G,v.)  =-<^V+  ^ 


^ip-iKa 

V 


2Tt 

ORDEP  £  1  (FIRST) 


(36) 


ay 


3^  2jac  \_724  -X>Svo 

—  6><?V o'y  +•  ^  ^  ^ 


(39) 


After  similar  substitution,  the  enthalpy  equation  in  the  wall  layer  yields 
the  following  hierarchy,  for 


4^  -fc^o 


ORDEP  1/fc 


-r 


4>,  = 


_  _v/ 

-  v°^ 


(40) 


(41) 
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Note  that  on  the  left  hand  side,  according  to  Eqs.  (28),  (33) 

'hv.j'Zi)  O 

(since  0  in  general).  Now,  according  to  Eq.  (36) ,  p0  -  p_(t) ,  60  that 

both  sides'  are  identically  zero,  as 

4i=4>.ct)=  -m  co/or-o 

Therefore,  Eq.  (41)  does  not  yield  any  new  information. 


ORDER 


C°  ( 2ER0TH) 

^  V.-r7SfoV,)~  -TTcfeVc, 


351 

fa  ^ 


ORDER  E1  (FIRST) 

M  U.+TTU4!,)  (r<t>  v,)  =  -TSftv.-y + r(<ftv,+c f,  v.) 

+  (T-0{ 

+  ■£-  /'-'at* +'§5>‘  ^ 

Pr  >»  "2JVJ  T  "E^p-  /  • 

(43) 

This  concludes  the  derivation  of  the  perturbed  equations  of  motion. 

The  results  of  lower-order  analysis  can  be  summarized  as  follows.  From  Eq. 

(28)  , 

foVo  -  =  F*  (*/*") 


while  from  Eq.  (2.87),  using  the  last  equation, 

vc  «  vc  ( x , t )  (44) 

Thus  Pa,  vD  and  F_  are  all  independent  of  y.  Further,  from  Eqs.  (31)  and 
(36)  clearly 

while  from  Eq.  (32) , 

Pi  *  Pi  (x,t)j  (45) 

so  that  both  pQ  and  pi  are  independent  of  y.  As  a  consequence  of  Eq.  (45), 

+ftho  =  4$  Cx-jf) 

One  may  now  proceed  to  solve  Eq.  (37)  directly  for  u_: 


(46) 


==>  « ^  ^  y 


(47) 
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with  the  boundary  condition,  (no-slip): 

ue  (x,o,t)  -  0  (4B) 

Of  course,  the  (x,t)  dependence  of  u0  still  remains  to  be  found.  However, 
its  dependence  upon  the  layer  coordinate,  y,  Is  found  to  be  linear;  this  result 
is  certainly  not  obvious,  and  has  several  important  implications. 

The  shear  stress  within  the  layer, 


is  obviously  nonzero  in  general,  while  being  independent  of  distance  from  the 
wall  at  any  given  station  (x,t).  Even  more  striking  is  the  vanishing  of  the 
viscous  dissipation  term  at  zeroth  order: 


which  leaves  in  the  zeroth-order  axial  momentum  equation  a  balance  of  inertial 
terms,  strictly,  cf.  Eg.  (38).  This  explains  physically  the  success  (up  to 
first  order)  of  modeling  this  family  of  injected  flows  by  assuming  rotational, 
inviscid  motions;  such  modeling  indeed  obtains  solutions  for  the  axial  velocity 
profile,  which  satisfy  the  no-slip  condition  at  the  wall  (r«l). 

It  further  appears  that  the  shear  stress,  Eq.  (49),  is  proportional  to  the 
first-order  axial  pressure  gradient,  while  being  inversely  proportional  to  the 
injected  mass  flux,  as  would  be  expected.  Of  eour  se,  SF| /^x  depends  on  FQ,  and 
one  expects  their  ratio  to  be  finite  at  the  limit  as  zero  injection  is 
approached. 

In  addition  to  the  foregoing  result  for  axial  velocity,  the  zeroth  order 
formulation  can  be  utilized  to  solve  for  the  y-dependence  of  the  other  dependent 
variables.  The  zeroth  order  continuity  equation  can  be  written  now  as 


Now,  since  FQ  and  are  independent  of  y,  hence,  one  may  split  the  foregoing 

equation, 


=  CoOst) 


ijl  (  _  1£l _ s-  (y. 

Vo  J  tjvj  (51 

where  C0(x,t)  is  a  common  separation  parmeter,  with  a  range  of  value-  uniquely 
corresponding  to  the  boundary  data. 


The  second  equation  yields 


l’za-v* 


where: 


Fi (y.x,t)  -  B0(x,t)  *  C_(x,t) .y4i  y 


B0(x,t)  m  Fx(0,x,t) 


g0(X,t)  ■  ^r  -  ^  . 

Similarly,  from  the  zeroth  order  radial  momentum  equation: 

tt/vt  ■t'FwVo  •=  c,  (*,*■) 
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w  ^  -  iw  (v'«F.+ 'af. )  =  -c,u,t) 


The  last  equation  can  be  integrated, 

V|C*J,)t,T)-B10‘«+  c‘~y°c?u+l  fe?  H&fat&Xp. 

*-*  J  2  fo  UK  77  y 

Bi-ViCOjX^)  . 

The  foregoing  results  for  F2  and  V2  yield  for  the  first-order  density: 

C&o  f©& iV^o  •+*  ^2Co  —  G< /W) 

-t-X.fe^rafc 

2  v.'  L  ax1  -tx  •?>:  J  “ 

Note  that:  ^(O.x.t)  -  (BQ  -  f’0B1)v0 

so  that  ^  is  uniquely  defined  and  an  additional  integration  constant  is 
necessary. 

Now,  following  the  definition  of  ^  and  since 
then: 

h^y.x.t)  -  -j5“^(y,x,t)  +  B2(x,t) 

Thus,  hjty.x.t)  is  second  order  in  y,  like  fjj  the  appropriate  boundary 
condition  is, 

B2(x,t)  «  h1(0,x,t)  -  h0(B0  "foBil/Fp 
The  axial  momentum,  after  dividing  through  by  Fc,  Eq.  (38)  leads  to: 

J-  _  X.  d-Cu*  .U*3£e_+3;-®l?-±r'&j.  J,® 

fi,  K>  ^  fi>  Tp.  *  V*- 

—  El.  DL^LLe  s-  _  | 
r=o  Try  '• 

Now,  according  to  the  foregoing  results: 

l  _  Co  , 


<3Vo  __  ^  _  I 

R  3^  -  1 


_ x !£  =  _x  |2i.y_& 

Fe  ”2^  'o£x.  F© 

=  -  p./ft’J 

Substitution  into  Eq.  (  41  )  yields,  after  some  manipulation: 

Tyi?"Sk.-i 


&u 
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where 


thus, 


(66) 


%  Mf-i  t  j  % 


and 


—  B*3Z°  *u 


(67) 


UjtO.x.t)  -  0 

satisfying  the  no-slip  condition  at  the  wall.  Thus,  the  perturbed  axial 
velocity  is  third  order  in  its  y-dependence,  and  the  corresponding  viscous 
dissipation  term  (unlike  its  zeroth-order  counterpart),  does  not  vanish  within 
the  layer. 

With  the  foregoing  derived  results,  the  zeroth  order  energy  equation  can 
be  shown  to  yield  merely  a  condition  conection  UQ  and  vc  timewise: 


U/v.  =  C4Ct) 


One  may  turn  now  to  the  first-order  energy  equation,  which  seems  to  yield  some 
simply  and  highly  useful  results  even  without  full  solution.  Equation  (43) 

0  written  in  terms  of  pressure,  reads: 

H  +  -r  1-  EM,  J  M  =  -r(p„v.y+  P,V»  P.v,l 

-Kr-Ou.'ig)  ■+ 

(69) 


Pr 


After  some  manipulation  one  obtains: 


i  <■' *°>rp°  + v^+v->  “-I  -"W  ^,oi 


The  first  bracketed  term,  after  using  Eq.  (68)  and  the  foregoing  results,  is 
simply 


V0  “  (Cl_vo<'0^ /^O 

The  enthalpy  term,  according  to  Eqs.  (58)  and  (59)  is: 


vaA  ^  ^  ;-vc 


(71) 


The  first  order  axial  velocity  gradient  car  be  found  from  Eq.  (67). 
Substitution  of  the  last  results  along  with  the  appropriate  expression  for  vlr 
into  Eq.  (70)  and  collection  of  equal  powers  of  y  yields: 
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{%  +7rR  O*  -  >v- 

<^t^V?c  i-  Ct||  }  V 

+{r?, Kfc*%5 -rft 

+  £?*  lKfc>  >-W3  =  0  • 

Compatibility  with  the  foregoing  derivation  (in  which  y  and  (x,t)  variable 
separation  was  implemented),  can  be  maintained,  provided  each  of  the  bracketed 
terms  in  Eq.  (72)  vanishes  identically.  The  resulting  four  compatibility 
relations(partial  differential) would  determine  the  behavior  of  the  wall  sublayer 
system  up  to  the  rirst  order  in  £  ,  the  small  perturbation  quantity.  However, 
a  total  of  three  undetermined  coefficients  should  arise  necessarily,  to 
accommodate  coupling  with  the  outer,  inviscid  (core)  flowfield,  through  inner/ 
outer  asymptotic  matching. 

THE  FIRST  ORDER  PRESSURE  PERTURBATION 

Of  particular  interest  in  the  present  analysis  is  the  pressure, 

P(y,x,t)  -  pQ(t)  +  £,Pi(x.t) 

which  is  a  directly  measurable  quantity.  From  the  axial  momentum  balance  in 
perturbed  form,  it  is  evident  that  the  rotational  ("inviscid")  coreflow  can  not 
sustain  a  first  order  term  like  '3p1/3x  herein;  the  lowest-order  axial  pressure 
gradient  effect  evolves  only  at  second  order,  or  £2pj-level.  This  is  clearly 
borne  out  in  the  analyses  of  Culiek,  and  others,  in  which  the  axial  pressure 
drop  is  proportional  to  YM02  (Mach  number  of  injection,  squared),  or  to  £,2 
according  to  the  convention  employed  here. 

This,  however,  is  not  what  is  observed  in  the  recent  injected  cold  flow 
experiments  of  Brown,  et  al  at  CSD/UTC;  the  measured  axial  pressure  profiles 
clearly  indicate  variation  of  order  MQ  _  £  ,  or  first  order. 

It  therefore  seems  that  the  viscous  wall  layer,  with  its  inherent  first- 
order  dissipative  processes,  impresses  this  axial  pressure  variation,  at  first 
order,  over  the  entire  cross  section  of  the  injected  channel.  Needless  to  say, 
this  would  involve  a  corresponding  variation  or  distortion  of  the  coreflow 
radial  and  axial  velocity  profiles  -  from  their  zeroth  order  representation. 

To  resolve  the  axial  variation  of  pi  by  the  wall  layer  formulation,  the 
second  compatibility  condition  in  Eg.  (72)  can  be  used,  corresponding  to  the  y-te 


N/.  j_  C-i  Co^e 

ax',  ~rr 


TPo 


(73) 


For  the  special  case  of  uniform  (zeroth  order)  injection  at  steady  state,  the 
presence  of  a  nonzero  first-order  pressure  perturbation  would  imply  physically  a 
corresponding  nonzero  perturbation  upon  the  mass  flux  injected,  i.e., 

Be  ■  F1(0,x,t)  /  0 

as  given  by  Eq.  (53).  Now  at  steady  state,  although  B „  is  expected  to  vary 
with  p1#  we  will  assume  for  simplicity  that  B0(x)  ■  b0  ■  const. 

With  the  foregoing  steady  state  assumptions,  Eq.  (73)  can  be  readily 
turned  into  an  ordinary  differential  equation,  for  0<x<L: 


d7rt 


o 


(74) 
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xhere  p^(x)  is  the  steady  state  first-order  pressure  perturbation;  the 
coefficients  are: 

l  ==  .  _  p;/^ 

Kw'vo  J  ^  TTfe  ns) 

Note  that  at  steady  state,  according  to  Eqs.  (50  )  and  (55),  respectively,  C0»F0 
and  Cj-FpV^  thus,  in  Eq.  (73),  C1-Covo*0. 

The  boundary  data  are, 

dpj/dr(0)»0,  and  p1(x»L)»pL  (76) 

The  solution  is  straightforward, 

dZ/dK  =  -\[b 7lb[  (77) 


tt  C^)  *=  f;  °  +■  I  Crtfykib?  *)  |  .  (78) 

This  concludes  the  derivation  of  the  injected,  viscous  wall  layer,  ud  to 
second  order.  Full  solutions,  namely,  matching  between  inner  and  outer 
expansions  will  not  be  attempted  herin.  Important  insights  are  obtained  already 
from  resolvinq  the  near-field  behavior  up  to  second  order,  in  terms  of  the  y- 
polynominals. 

DISCUSSION  OF  RESULTS 

To  facilitate  comparison  with  the  experimental  data  reported  by  Brown,  et 
al,  one  may  form  the  normalized  axial  pressure  differential, 

^ = at =spj p(K=o)-  *  i 


(79) 

This  axial  pressure  differential  expression  is  used  to  correlate  the 
experimental  data  of  Brown,  et  al,  as  shown  in  Fig.  3.  Clearly,  the  measured 
pressure  orofile  is  correlated  very  well  by  Eq.  (79),  which  is  obviously 
superior  to  the  expression  attributed  to  Culick,  shown  as  well. 

It  should  be  pointed  out  that  a  single  point  of  the  data  (x:  Ap^)  has  been 
utilized  to  obtain  a  scale  for  the  comparison  (this  is  necessary,  since  no 
physical  input  is  available  regarding  the  value  of  Bp,  the  perturbed  injected 
mass  flux,  necessary  for  defining  b^,  bD) ,  along  with  p0«F_  ■  v_  ■  1,  and  If' 

1.4.  Suppose  now  that  Km*l,  and  we  select  a  value  of  Bo“60.  (This  is  based  on 
some  trial  and  error  -  but  shows  how  the  correlation  was  obtained  without  any 
regression  analysis);  then, 

bo«l/Bo»l/60,  bj»l/  1  Bo«l/1.4x60 ,  •  0.012 

m  Vf*  Bo  *  °-014 

Two  important  observations  are  therfore  demonstrated:  (1)  axial  pressure 
variation  to  lowest  order  is  0(£.),  and  is  governed  by  the  dissipative  wall 
layer  processes,  as  shown  in  the  rigorous  analysis  herein.  The  behavior 
obtained  in  x  differs  from  than  the  parabolic  pressure  drop  formula  of  Culick 
(1),  and  (2)  one  need  not  invoke  local  turbulence  generation  or  turbulence 
encroachment  upon  the  surface  to  explain  the  departure  of  measured  Pj^  from  a 
laminar  behavior. 

Another  oropertv  of  interest  is  the  wall  friction  coefficient,  or 
dimensionless  wall  shear  stress,  # 

Cf.  =  ~c2/W  u*  -  -K*  / kfu? 
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here  u*  denotes  the  mean  axial  coreflow  velocity.  Using  dimensionless 
invention  employed  herein,  along  with  the  wall  layer  coordinate, 

s  u»2x  was  used,  for  a  cylindrical  port,  and  subscript  zero  denotes  zeroth 
rder  convention. 

Now,  from  Eqs.  (47)  and  (77), 

° v  (8d 

here  the  first  square  root  term  is  of  order  unity.  This  parameter  is  plotted 
gainst  l/2x  (which  denotes  the  ratio  of  blowing  to  mean  axial  velocity)  in  Fig. 
.  A  nearly  linear  relationship  is  obtained,  using  the  foregoing  coefficient 
alues.  In  comparison,  the  data  obtained  by  Olson  and  Eckert  [6]  is  considered, 
ief.  6  includes  a  plot  of  the  ratio  of  (axial  pressure  gradient)/(meen  dynamic 
ixial  head)  vs  vQ*/u0*  ■  l/2x.  This  obtains  an  almost  linear  correlation,  as 
rould  be  exoected  from  a  parabolic  pressure  drop.  The  slight  curvature  however, 
larticularly  apparent  at  small  values  of  l/2x  <  0.01,  can  be  followed  only  with 
he  present  formulation,  not  with  any  parabolic  pressure  profile.  Thus,  the 
first  order  pressure  distribution,  obtained  from  he  viscous  wall  layer  analysis, 
iarees  well  with  the  measured  data  of  Brown,  et  al  18],  while  the  associated 
*all  friction  coefficient  follows  the  same  trend  as  that  measured  by  Olson  and 
ickert  161. 


CONCLUSIONS 

A  derivation  of  the  viscous  wall  layer  regime  has  been  presented, 
pertaining  to  injected  flow  in  an  axial  porous  tube,  in  simulation  of  interior 
solid  propellant  rocket  flows. 

Solutions  for  the  radial  coordinate  (or  y-dependence)  of  all  the  deoendent 
variables  u?  to  te  second  order  have  been  generated,  in  polynominal  form.  The 
(x,t) -dependence  is  defined  in  terms  of  a  relatively  simple  partial  differential 
system. 

Particular  results  of  the  analysis  for  the  SDecial  case  of  steady  state, 
sre:  (1)  the  first  order  pressure  perturbation  was  solved  for  and  its  axial 
distribution  is  given  explicitly;  this  term  is  entirely  due  to  the  laminar 
dissipative  wall-layer  processes,  and  (2)  the  blown  wall  friction  coefficient 
was  likewise  defined.  Both  results  correlate  well  the  available  experimental 
data.  Finally,  (3)  the  zeroth  order  axial  velocity  distribution  within  the 
layer  is  linear  radially;  thus,  to  lowest  order,  viscous  dissipation  is 
negligible  in  the  axial  momentum  balance.  This  indicates  why  inviscid, 
rotational  solutions  (such  as  those  of  Culick  (13  and  others,  chosen  so  as  to 
Batisfy  the  no-slip  condition  at  the  wall)  are  so  successful  in  representing 
these  family  of  flows  -  up  to  first  order. 

NOMENCLATURE 

Ht,  A  »  nozzle  throat  area  and  port  exit  area,  respectively 

a  *  adiabatic  velocity  of  sound 

Cj  ■  wall  friction  coefficient,  Eq.  (80) 

cw,  Cp  •  isochoric  and  isobaric  specific  heats  (J/kg-K) 

F  »  radial  mass  flux  (dimensionless) 

G  •  axial  mass  flux  (dimensionless) 
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n 

L 

W 

P 


r 


S1.2,3 


«  thermal  enthalpy,  dimensionless 

■  ratio  of  inverse  Reynolds  number  and  Mach  number  squared,  Eq.  (21) 

■  chamber  length 

*  Mach  number 
«  pressure 

■  Prandtl  number,  Eq.  (7) 

■  channel  radius 

■  injected  Reynolds  number,  Eq.  (7) 

•  radial  coordinate 

■  "source’-terms  in  the  eauations  of  motion  for  coreflow,  Eqs.  (12)— 
(14) 


SR0  ■  Strouhal  number,  injected,  (Eq.  24) 

t  “  time  (dimensionless) 

U0  «  parameter  defining  (x,t)  -  variation  of  wall  layer  axial  velocity 

component 

u,  u  *  axial  velocity,  and  mean  axial  coreflow  velocity  respectively 

v  *  radial  velocity  component 

x  «  axial  distance  (dimensionless) 

y  ■  radial,  magnified  wall  layer  coordinate,  oerpendicular  to  surface, 

Eq.  (18) 

Greek  Symbols: 

■  Cp/Cv  specific  heat  ratio 
«  difference,  increment,  Eq.  (79) 

»  length  scales,  Eqs.  (22) — (24) 

■  small  perturbation  quantity,  Eq.  (19) 

■  thermal  conductivity  of  gas  (air),  J/K-m-s 
«  viscosity  coefficient,  kg/m-s 

■  density 
Subscripts,  Superscripts: 

(  )D  ■  denotes  zeroth  order  (perturbation) 

(  )j  »  denotes  first  order  perturbation 

(  )*  ■  denotes  dimensional  quantitv 
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WITHIN  THE  VISCOUS  SUBLAYER  (STONES  LAYER)  WHERE  PRIMARY  COMBUSTION  OCCURS: 


PRECUENCY-DEPENDENT  SURFACE  DC-COMPONENT:  ACOUSTIC  STREAMING, 

HEAT  TRANSFER  (P1RST  ORDER),  WHEN  U,  V  NOT  OUT  0F  PHASE,  NET 

with  phase  relative  to  x-momentuh  transpep  occurs. 

ATTENDANT  PERTURBATION.  (SECOND  ORDER), 


nonsteady  CONDENSED 
PHASE  (THERMA. 
RELAXATION,  DynA-IC 
BURn) 


Figure  1.  Schematic  of  the  inner  motor  nonsteady  combustion 
process  simulated  by  the  injected  cold  flow  analysis. 


WITH  INJECTION 


Figure  2.  Schematic  of  the  injected  wall  layer  region 
showing  the  magnified  coordinate. 
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Fioure  4.  Kail  friction  coefficient 
vs  reduced  blowing  rate  showing 
slight  curvature  at  small  blowing 
fraction.  This  agrees  well  with  the 
results  of  Olson  and  Eckert  16),  and 
could  not  be  explained  bv  the 
inviscid  solutions  of  Culick  11). 
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